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DIGITAL COMPUTER PROGRAM FOR ANALYSIS OF CHUGGING INSTABILITIES 

by John R. Szuch 
Lewis Research Center 


SUMMARY 

A digital program, designed to generate chugging stability limits for bipropellant 
rocket engines, has been developed and is described herein. The program is written 
in FORTRAN IV language for use with the IBM 7094 computer. The program computes 
combustion delays, gas residence time, characteristic velocity, and other steady- state 
parameters required for solution of the characteristic equation. The characteristic 
equation is solved by the computer for critical values of the injector pressure drops and 
the chugging frequency. Stability limits have been generated for a specified flox-methane 
engine system. Results have been obtained for two injector configurations. For each in- 
jector configuration, limits for throttling at constant mixture ratio and the effects of 
varying the fuel injector area on stability were determined. Digital results were com- 
pared with available analog data to test the program. For both injector configurations, 
the predicted throttling limit was about 5 to 1 because of higher frequency (500 Hz) insta- 
bilities, characteristic of the chugging model. Typical computer listings and printouts 
of results are presented. 


INTRODUCTION 

One of the problems often encountered in the development of liquid-propellant rocket 
engines is the occurrence of low-frequency instabilities, commonly referred to as chug- 
ging. Chugging is caused by a coupling of the propellant feed system with the combustion 
dynamics in such a way as to reinforce any disturbance in pressure or propellant flow. 

Many analyses of chugging have shown that the controlling factors in determining 
low-frequency stability are the time delays associated with the atomization, vaporization, 
and combustion processes. Without these delays, the chamber- injector-feed system 
would be inherently stable. 

The validity of a double-dead-time model (ref. 1) in analyzing chugging has been 
demonstrated in references 2 and 3. Previous analyses involved the use of a single- 


dead-time model (ref. 4) with combustion delay, calculated for the propellant having 
the longer drop lifetime, applied to both propellants. The double-dead-time model re- 
quires the application to each propellant of its respective time delay and results in sig- 
nificantly different stability characteristics. Experimental data, reported in references 
2, 3, and 5 to 7 have been matched using the double-dead-time model. 

The implementation of chugging models has traditionally been done on the analog 
computer. However, as the models became more sophisticated and as engine configur- 
ations and propellant combinations changed, the need for a high-speed digital stability 
program became more evident. 

A digital program, designed to generate stability limits for bipropellant rocket en- 
gines, has been developed and is described herein. The program was written in FOR- 
TRAN IV language for use with the IBM 7094 computer. Linearized feed-system im- 
pedances are evaluated by the program at each frequency of interest. In addition to 
chamber and injector flow dynamics, the program can handle combinations of pumps, 
valves, manifolds, and lines. For any selected engine configuration, the program com- 
putes the required steady-state engine parameters for solution of the characteristic 
equation. Available vaporization (ref. 8) and drop- size (ref. 9) correlations are used 
to calculate combustion delay times. The characteristic equation is solved by the com- 
puter for critical values of the injector pressure drops and the chugging frequency. 

Stability limits have been generated for a specific flox (82. 6 percent fluorine) - 
methane engine system. An expander cycle with a turbine -bypass throttle has been 
assumed. Results have been obtained for two candidate injector types. For each injec- 
tor configuration, limits for throttling at constant mixture ratio and the effects of vary- 
ing the fuel injection area on stability were determined using the digital program. Digi- 
tal results were compared with available analog data to test the program and the lineari- 
zation procedures. 

Although the existing computer program is designed for the case of a completely 
vaporized fuel, the program is easily modified to handle liquid- liquid propellant com- 
binations. The steps required for this modification, together with other aids to the user, 
are included in an appendix to this report. 


ANALYSIS 

Derivation of Characteristic Equation 

The combustion process in a bipropellant rocket engine is still not clearly under- 
stood, although what takes place is qualitatively known. Figure 1 shows, schematically, 
the assumed model for a bipropellant rocket combustor. Initially, the propellants are 
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nozzle 

Figure 1. - Schematic of bipropellant combustor. 


injected into the combustion chamber. The propellants then atomize, vaporize, and 
mix; they then react to produce hot gases. Their flow rates are determined by the con- 
ditions upstream and downstream of the injector element and the injector geometry. 
These processes are gradual and continuous. For mathematical expediency, however, 
they are treated in a discontinuous manner. A time interval between injection and a sud- 
den conversion to vaporized propellant is assumed and referred to as the vaporization 
time delay cr y . The time interval between the conversion to vaporized propellant and 
the conversion to burned gas is similarly referred to as the gas-phase mixing delay <7 . 
For this analysis, the time delays are assumed to be invariant with time. With these 
assumptions, the following equations relate the rate of change of burned products to the 
injected propellant flow rates: 


^obW = w iQ (t - a yo - aj 

Wfb(t) = w if (t - d vf - a m ) 

(All symbols are defined in appendix A. ) 

If the burned products behave as a perfect gas, 



( 1 ) 

( 2 ) 


( 3 ) 


Assuming a small disturbance in the rate of change of burned products Aw^ and using 
the conservation of mass yields 
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(4) 


d_ 

dt 


[w c (t)] » W ob (t) + Wfb(t) + Aw b (t) - w n (t) 


For a choked nozzle, 


V‘> * 


PqWVc 

c*(t) 


(5) 


Equations (1) to (5) can be combined to give the following equation relating total pressure 
at the nozzle throat to the injected propellant flow rates and the disturbance: 


-X. c - c *(t) £ 

A t g c dt 


p c (t) 


R T it) 
gc c v f 


+ p c w - 


— [wio(t - °To> + "if* 4 - a Tf } + AW b (t) ] 


A t g c L 


( 6 ) 


where 


a To = °Vo + a m 
a Tf = °Vf + °m 


Equation (6) is nonlinear and must be linearized to obtain the desired analytical solution. 
Equation (6) was linearized by assuming small perturbations in the system variables 
about a steady- state operating point and neglecting products of perturbations: 


c* 


dP 


c(t) 


\e c (R g T > c «t 


+ ? c(t) Pio (t - a To> + "'ll* 4 

A t g c 


(J Tf ) + Aw 


v b (t)] 


' w o + w f) 


A tS, 


c*(t) 


t s c 


A t g c (RgT)^ dt 


(R g T) c (t) 


(7) 


The coefficient c*/[ A t g c (RgT) c ] is the gas-residence time and is invariant with 
time. It is generally accepted that for small perturbations from steady state, the cham- 
ber temperature T c remains nearly constant. Therefore, the last term on the right 
hand side of equation (7) was set equal to zero. Equation (7) is then rewritten as 
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dP 


_£(fi + 5 c 

dt c 


(t) = A^~ F io(t ' CTto) + ^ if(t ‘ CTTf) + A ^ b(t) J 


' w o + w f > 


A t^c 


c*(t) (8) 


The theoretical characteristic velocity c^ is a known function of the chamber mix- 
ture ratio; that is, 



(9) 


Om. W io (t ~fTo) _ *ob (t) 
F w if (t - a Tf ) Wfj^t) 


( 10 ) 


Linearization of equations (9) and (10) yields the following equation relating perturbations 
in c* h to perturbations in the injected propellant flow rates: 



Assuming that the actual characteristic velocity c*(t) is related to the theoretical value 
by a time- invariant efficiency, 


c*{t) = T 7 C * c* h (t) 


( 12 ) 


equations (8), (11), and (12) can be combined to yield 


_ dP (t) ^ _ 

C . T» _ -V.:. / 1- ~ \ . iP.:» It- rr \ . C Aw^(t) 


e 


g 


dt 


+ P c^ = XW io^ - "To* + F ^if^ - CT Tf) + 


A t^c 


(13) 
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where 



The Laplace transform of equation (13) is 


(ejs- 


1)P C (S) = Xw io (S)e 


' a To S 


Fw if 


-a„ f S 

(S)e Tf +. 


c* 


A t§c 


Aw, (S) 


( 14 ) 


(15) 


(16) 


In general, the perturbations in the injected propellant flow rates can be related to 
the perturbations in the chamber pressure at the injector face P . by 


w i0 (S) 


Z so< S > 


P ci (S) 


(17) 


w [f (S) - - 


Z sl (S) 


P ci< S > 


(18) 


where Z gQ (S) and Z g ^(S) are the linearized output impedances of the oxidizer and fuel 
feed systems, respectively. The impedances are evaluated at the injector face and in- 
clude the impedances of the injector, valves, pumps, etc. 

The chamber pressure at the injector face P ci and the chamber pressure at the 
throat P £ can be related by considering the time rate of change of momentum of the 
flow between the injector and the nozzle inlet. By equating the pressure force and in- 
ertia force on the gas, one obtains 


w. 


Pci Pc 


A g 
c & c 


< V C - v ci ) = 


v p 
crc 


v ci p ci 


Sr 


(19) 
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where p ci and p c are the static chamber pressures at the injector face and nozzle in- 
let, respectively. Equation (19) may be written in terms of the Mach number such that 


Pci _ 1 + y c M c 

P c 1 + r c ^ci 


( 20 ) 


By definition, the total pressures P c and P ci are related to the static pressures p c 
and p c p respectively, by 


P c = Pc 


1+1 



y c /{y c -V 


( 21 ) 


Cl 


Cl 



-I y c /( r c' 1) 


(22) 


Combining equations (20), (21), and (22) and assuming that is small (large 

chamber) results in 


P ci 


1 + r c M c 



-,r c /(r c -V 


= K 


(23) 


By assuming that the pressure ratio K c is time- invariant and that there is no total 
pressure loss in the converging section of the nozzle, equations (16), (17), (18), and (23) 
may be combined to yield the following transfer function: 


c* 

p c < s > _ ^ _ 

A ^ (S > Vs ♦ 1 ♦ ?T0 S + ,- _C Tf S 
g z so‘ s > V s * 


(24) 
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The transfer function describes the response of the total pressure at the nozzle inlet P c 
to the disturbance in the rate of change in burned products Aw^. 

The characteristic equation, describing the stability of the bipropellant system, is 
obtained by setting the denominator of equation (24) equal to zero. The result is 


XK c e - CT To S + e "°Tf S 

1 _fso®__ + V S L 


(25) 


It can be shown that equation (25) applies regardless of where the disturbance is intro- 
duced into the system. 

Figure 2 is a block diagram of the linearized bipropellant system described by equa- 


Oxidizer 
injector 
and feed- 
system 
dynamics 



Fuel in- 
jector and 
feed -system 
dynamics 


Figure 2. - Block diagram of chugging model for general bipropellant engine system. 


tion (25). If all of the roots of equation (25) have negative real parts (S = X + jco, X < 0) 
the system is stable. The first appearance of a root on the imaginary axis (S = jw) de- 
fines the boundary of stability. 


Solution of Characteristic Equation 

By separating the impedances Z gQ and Z g j into their respective real and imagi- 
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nary parts and letting S = jo>, the characteristic equation (eq. (25)) can be separated to 
give 


! XK c ^o cos ( u,or q' 0 ) " XK C >o sin(u)arp o ) FK c # f cos(wa Tf ) - FK C > f sin(cocr Tf ) 


gt 2 + j ^ 
* o + •'o 




(26) 


XK C cos(cua To ) + XK C dt Q sin(coa To ) FK c > f cos(cucr Tf ) + FK c # f sin(u>cr Tf ) 
“ v <*>; + — ~ 


(27) 


where 


St, 


+ 1 


•'o 


= Z so^> 


&f + iS f = Z sf (jco) 


Specifying angular frequency o> allows equations (26) and (27) to be solved for two 
selected critical parameters. The real parts gt Q and were selected because they 
could be related to the injector flow resistances d AP/dw which are independent of fre- 
quency. Specification of a parameter other than frequency requires the iterative solution 
of transendental equations involving frequency. 

Appendix B outlines the steps required to solve equations (26) and (27). The pro- 
cedure followed results in the following quadratic equations relating the critical value 
of the real part of the oxidizer impedance &' 0 (u') to the imaginary parts of the feed 
system impedances S Q (o)') and > f(cu’) (the engine gains, delay times, imaginary parts, 
etc. , are evaluated at the selected operating point): 


3t\ 


2XS f 

- w’ovpf) - — — — 

F 


To w Tr 


Y - = ( 1 + - ,20 |)j +*o{ -XK c 

jcos(co r a To ) - w'$ g sin(u)’a To )j| + - > Q XK c cos(w’a To - a Tf ) 


-^-22 


[ Sin(a, ’ a To) + w '®g cos ^'*To>] ' ^ i 1 + “' 2Q l) 


• X 2 ^ 


= 0 


(28) 
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where 


Y = c o’0 cos(co’a Tf ) + sin^'o^) 

The solution of equation (28) is used to compute the corresponding critical value of 
the real part of the fuel impedance dip The equation, which is derived in appendix B, 
is 

_1_ = VV> 0 -XK c sin ^So ) _ 

FK c ^q sin(w T CT Tf ) - F J Q K c cos(co'a Tf ) - XK^ cos(a;'a To ) + co'0 g 7 o > f - 7 f ^ 

(29) 

The values of 3?^ and dtp which satisfy equations (28) and (29) at frequency co’, 
must be related to the physical injector and feed systems being studied. In addition, the 
imaginary parts s and J ^ must be computed for the particular feed system and 
operating point at each frequency before solving equations (28) and (29). The discussion 
of how this is done is contained in the next section. 


Evaluation of System Parameters 

The primary function of the computer program is to test the stability of a specified 
engine configuration operating at a selected thrust level or chamber pressure P . The 
program may be used to study the effects on stability of changing the engine geometry at 
fixed thrust, throttling the engine with fixed geometry, or combinations of both. Regard- 
less of the type of study, solutions of equations (28) and (29) to find the critical feed sys- 
tem impedances dt' Q and require prior calculations of the system parameters at 
the specified operating point. These calculations are carried out by the computer pro- 
gram. 

Engine gains . - The oxidizer and fuel gains X and F, respectively, are computed 
using equations (14) and (15). A curve of the theoretical c* against mixture ratio for 
the specified propellant combination is required. For this study, equilibrium properties 
are used with the propellant temperatures equal to full thrust values. Although a dynamic 
sensitivity of c*^ to only mixture ratio is assumed, a static sensitivity to both mixture 
ratio and chamber pressure is used to determine the steady- state value of c*^. The 
digital program computes both c^ and 3c*^/ 3(0/F) for known values of P c and O/F. 
Required input information also includes the throat area A^. and combustion efficiency 
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V c *• The parameter r ) c , * is computed at each operating condition and a discussion of 
the calculation follows. 

The chamber momentum loss parameter K c is calculated using equation (23). The 
specific heat ratio y is determined from specified values of P and O/F. The Mach 
number at the nozzle entrance M c is calculated from 


M c=- 


(y c +l)/2(y c -l) 


e cYc + 1 J 


(30) 


where e is the nozzle contraction ratio or A /A.. 

C C L 

Co mbustion efficiency . - The combustion efficiency 77 c * is a measure of the degree 
of mixing and vaporization of the propellants within the chamber and nozzle. Work de- 
scribed in reference 8 was directed at relating the efficiency due to incomplete vaporiza- 
tion T? va p to the chamber geometry and operating conditions. For that study, perfect 
mixing of propellants was assumed. It is here assumed that 


^c* ^vap^mix 


(31) 


The following equations were derived in reference 8 to compute rj. 


vap’ 


_ *[£(§)] <*„***, 


'vap 


(32) 


w Q + w f 


where 

C fraction of oxidizer mass vaporized within chamber nozzle; function of l 

gen , o 

& fraction of fuel mass vaporized within chamber nozzle; function of l , 

gen, f 

^ functional relation of c£ h and mixture ratio 

and where Z is a generalized length parameter 
& on 
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0.66 


■J, 


0.83 Z 


l 


gen, o 


_0. 44 


0. 83 l 


' gen, f 


0. 44 


0. 22^0. 33 1 
e c y 


^2. 07x10 


0. 22„0. 33 1 
€ c y 


*. n 


1.45," v 0. 75 


cr, o/ 


V7. 62X10" 


P-) 7— °j 

^30. 5 y \100j 


0.35, 


H 


{3. 26x10" 


\ 0 . 66 


l 2. 07x10 


» t > \0. 4/ 


1.4 5/ - 


1 -■ 


cr, fy 


V7. 62X10 


-5y 


30. 5, 


0. 75 Sjf \ 0- 35/ 


100 


H f 


0.8 


26x10 J 


(33) 


0.8 


(34) 


Calculation of the generalized length parameter for each propellant requires knowl- 
edge of the propellant drop radius r. For liquid-liquid systems, the drop radii are 
usually assumed to be only functions of the injection element sizes (ref. 8). Experi- 
mental data presented in reference 9 for concentric tube injectors with gaseous fuel in- 
jected through the annuli indicated that the mean oxidizer drop radius r Q is also pro- 
portional to the square-root of the injection momentum ratio; that is, 


rw v ^' ® 

WJ^ 01 


(35) 


(concentric-tube injector, gaseous-fuel). 

Using experimental data presented in references 2, 3, 5, and 7, a proportionality 
constant K r of 0. 118 was calculated for liquid-oxygen droplets. However, references 
8 and 10 indicate the following effect of propellant properties on drop size: 

/ * x 0.25 

r cc j (36) 


where 

p, viscosity of liquid 

6 surface tension of liquid 

p density of liquid 
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To compute the drop radius for propellants other than liquid oxygen, the experimen- 
tally determined proportionality constant (K r = 0. 118) must be modified to include the ef- 
fect described by equation (36). The resultant expression for the drop radius of any 
liquid oxidizer droplet, injected through a concentric element with gaseous fuel injected 
through the annulus, is 


r o = 0- 118 d 


VA 2 
V 2 ^ 0 , 2 Pa l 


°- 25 / 5 „ v \ 0 ' 5 

o o 

Vf, 


(37) 


For this study, concentric tube injection with a gaseous fuel is assumed with the 
drop radius from equation (37) used in equation (35) to calculate the generalized length 
for the oxidizer. For the gaseous fuel case, is equal to 1. 0 with the efficiency due to 
incomplete vaporization of the oxidizer computed from 


^vap 



(gaseous fuel) 


(38) 


where denotes the functional relation of c* h to mixture ratio. 

Since the introduction of a swirl to the injected oxidizer stream will tend to enhance 
the atomization of the oxidizer, equation (37) was multiplied by 0. 448. This factor was 
based on available unpublished data for injectors having swirl elements. 

A suitable mixing model for calculation of v m i x was n °t available. Experimental 
data would be required to allow the use of a mixing model such as the one proposed in 
reference 11. For this study, equations (38) and (31) were used to calculate a value of 
77 m ix from an assumed T7 C * and a calculated 77 vap at full thrust conditions. This value 
was assumed to be constant over the entire range of studies performed. 

Comb ustion de lay times. - The chamber gas residence time was computed from 
specified chamber and nozzle geometry and combustion properties as follows: 


^c* c th / V \ 

g Vc \v/ c 


(39) 


Many analyses of chugging (refs. 12 and 13) have shown that the most important fac- 
tors in determining low-frequency stability are the time delays associated with the atomi- 
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zation, vaporization, and burning of the propellants. Experiments described in refer- 
ence 13 were aimed at defining the static relation between dead time and variables such 
as chamber pressure, mixture ratio, and injection velocity. The work described in ref- 
erences 14 and 15 was based on a dynamic sensitivity of combustion delay to chamber 
pressure and injection velocity, respectively. This study will consider only a static sen- 
sitivity of combustion delay times to the operating conditions. 

Studies reported in reference 16 equated the vaporization delay time a v to the time 
required to vaporize 50 percent of the mass of an injected droplet. Time histories of in- 
jected droplets (presented in ref. 8) indicate that the average droplet velocity over this 
length is approximately equal to the injection velocity. Using this information enables 
one to calculate the vaporization time delay from: 


a 


v 



(40) 


where v. is the injection velocity of the propellant droplet. 

The length required to vaporize 50 percent 1 5Q is calculated by solving equa- 
tion (34) for the chamber length l with the generalized length parameter equal to 0. 0699 
meter (ref. 8). This value corresponds to a fraction vaporized C of 0. 5 and applies 
for all propellants. The nozzle term in equation (34) is neglected for this calculation. 

The resulting expression for l is 


0. 44 


l 5Q = 0. 0699 



1. 45 


30. 5 


0. 75 


'_ufY 

100 / 


0. 35/ 


H 


0.8 


v3. 26x10" 


0.66 


2. 07x10' 


(41) 

For the case studied, the vaporization time delay was computed for the liquid oxi- 
dizer only. The gaseous fuel vaporization time was set equal to zero, since the 
fuel is assumed to be completely vaporized before injection. 

The gas-phase mixing and reaction delay cr m was related to the distance from the 
burning zone 7 5 q to the nozzle throat and the gas velocity at the nozzle inlet. Based on 
stability data from references 2, 3, 5, and 7, the curve shown in figure 3 was generated 
for 


cr = y\ 
m 1 


10 < ; c + *n- T 50> 


sec 


(42) 
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where the gas velocity v c was computed from 

v c = M c Vg c (^R g T)c m/sec (43) 

and S? represents the functional relation. 

Fe ed syste m i mpedances . - Solution of equations (28) and (29) requires evaluation 
of the imaginary parts of the feed system impedances as functions of frequency. In addi- 
tion, the critical values of &' o and 9t'^ must be compared with the operating point val- 
ues to determine stability. Therefore, the computer program must, for a specified feed 
system configuration, compute the real and imaginary parts of the feed system impe- 
dances at the specified frequency co’. This calculation, in general, will involve (1) the 
breaking up of the feed system into elements, each having a flow impedance Z^ with 
real and imaginary parts cy. and /3-, respectively, (2) the manipulation of these ele- 
ments (if necessary) into series and parallel combinations, and (3) the stepwise reduc- 
tion of these combinations using conventional reduction techniques (ref. 17) to obtain the 
real and imaginary parts of the feed- system impedance. 

The resulting values for the series combination of impedances Z^ and Zj (i. e. , 
z k = z i + z j) are 


cv k = Oj + <Vj (44) 

% = /3j + jSj (45) 
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For the parallel combination of Z. and Z- = Z-Zj/(Z. + Zj)J the real and imag- 
inary parts of Z^ are 


a k = 


_ aj (a? + e?) * or, (of + pf) 


+ a/^r + + #j)‘ 


(46) 


_p]( a i +f] i) +i 3 i(°'! +f3 !) 

P k 2 “ 2~ 

(c^ + + /3.)^ 


(47) 


Equations (44) to (47) form the basis for calculating the feed-system impedances at 
the frequency of interest regardless of the feed-system configuration. 

Using a lumped-parameter approach allows the liquid feed systems to be divided into 
the following elements: 

(1) Lossless lines having only momentum pressure drop. The flow impedance con- 
sists of a = 0 and 0 = col /Ag c , where Z is the line length and A is the cross-sectional 
area of the line. The electrical analog of this element is an inductance. 

(2) Orifices and valves having no momentum pressure drop and flow impedance con- 
sisting of a = 2 AP/w and 0 = 0, where AP and w are the steady- state pressure 
drops across the element and flow, respectively. The electrical analog of this element 
is a resistance. 

(3) Storage volumes having a shunt impedance consisting of a = 0 and 

/3 = -Bg /u>pV , where B is the liquid bulk modulus and V is the storage volume. The 
c g 

electrical analog of this element is a shunt capacitance. 

(4) A pump having no momentum pressure drop or storage. The pump is described 
by a normalized pressure-rise characteristic; that is, 


■= 

iP P ■= K ip* P N p 


(P 


P 


K 2p w p 

N 


> 



(48) 


For the frequencies of interest, the pump speed is assumed to be invariant with 
time, and the pump impedance consists of a - - K i p K 2 pNp( dip p / dq > p ) and 0 = 0, where 
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Np is the steady-state pump speed and dtp^/dcp^ is the slope of the pump characteristic 
at the operating point. The negative sign in the equation for a arises from the fact 
that the AP in equation (48) is a pressure rise, and the impedance d AP/dw is defined 
in terms of a pressure drop. It should be noted that operation is usually limited to that 
portion of the characteristic having a negative slope, giving an electrical analog of a 
positive resistance. 

Figure 4 is a schematic of the expander cycle engine to be studied. In this cycle, 
the fuel passes through a cooling jacket, which surrounds the combustion chamber. The 



Figure 4. - Schematic of bipropellant engine system using expander cycle. 


fuel is vaporized in the jacket and is used to drive a turbine before being injected into 
the combustion chamber. For the system shown, geared pumps are driven by the tur- 
bine supplying high-pressure propellants to the feed system. A turbine-bypass valve 
is shown; its function is to adjust the engine thrust by varying the amount of fuel 
passing through the turbine (hence, speed). A control valve is located in the liquid oxi- 
dizer system to maintain the desired mixture ratio in the chamber for all thrust levels. 

To test the computer program it was decided to study two oxidizer feed systems 
that would represent two degrees of difficulty in the evaluation of the oxidizer feed sys- 
tem impedance. Figure 5 shows an impedance representation of the two systems con- 
sidered for this study. Figure 5(a) shows a conventional single- orifice injector supplied 
by a manifold volume (capacitance). The injector elements consist of several long in- 
jection tubes (inductances) with pressure drop (resistance) controlled by orifices im- 
mediately upstream of the injection tubes. The tubes are long enough to guarantee 
steady entrance conditions at the chamber entrance. With a concentric-tube injector, 
the injector tubes are surrounded by annuli through which the gaseous fuel is injected. 
The oxidizer manifold is fed by a control valve (resistance) whose area is adjusted to 
maintain the mixture ratio. A pump discharge volume (capacitance) is located between 
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Figure 5. - Impedance representation of two oxidizer feed systems. 


the pump (resistance) and the control valve. The pump is assumed to be fed by a tank 
(battery) and a lossless suction line (inductance). 

Figure 5(b) shows a dual- orifice oxidizer injector and feed system. The conven- 
tional injection tubes and orifices are fed by two flow components. After discharging 
from the pump into a pump discharge volume, the flow is split. The bulk of the flow 
passes through a mixture-ratio control valve and into a secondary-flow manifold (capac- 
itance) which in turn feeds several secondary (axial) orifices. The remaining portion of 
the flow is collected in a primary manifold (The capacitance is lumped with the pump 
discharge capacitance. ) and distributed to a number of tangential orifices. For this 
study, each axial orifice had three associated tangential orifices. The tangential flow 
mixes with the axial flow in an annular chamber imparting a swirl to the injected flow 
stream. The vortexing action of the tangential flow (ref. 18) maintains the swirl at ex- 
tremely low thrust (flow) levels enhancing the mixing of the oxidizer stream and the 
gaseous fuel. For a particular thrust level, the area of the axial orifices is adjusted 
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statically in the computer program to satisfy the following assumed vortex character- 
istic: 


w 


sec 


— 

w o 



1. 1259 


( 49 ) 


Appendix C outlines the evaluation of the real and imaginary parts of the oxidizer 
feed systems shown in figure 5. At a particular operating point and frequency, the im- 
pedance Z gQ can be expressed as 


Z so^) = ^o + j-'o 


# Q = a + s/(cjo') > 


= u)'b + oo') 


(50) 


where 

a oxidizer injection orifice resistance 2 AP/\fr 

b injection tube inductance 

The imaginary part is used in the solution of equations (23) and (24). The real 
part £ Q will be compared with the critical value 3t' Q resulting from the solution of 
equation (23). By assuming that the injector geometry allows changes in the orifice 
resistance a without affecting either the parameters required to solve equation (28) 
or the rest of the feed system impedance the critical value of a can be found 

from 


a’ = dt' Q - j/(u>’) (51) 

where j/(cc’) is calculated at the selected operating point as outlined in appendix C. 

Although the same procedures hold for the evaluation of the fuel impedance Z g j, 
some assumptions are required to conveniently handle the system shown in figure 4. It 
was physically realistic and mathematically expedient to assume a choked turbine and 
bypass valve. Assuming steady conditions upstream of the turbine, the fuel injector 
can be considered to be fed by a constant fuel flow. Since the fuel is vaporized in the 
cooling jacket, the compressibility of the fuel must be considered in the evaluation of 
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Z gf . The following equation can be used to compute the fuel flow rate through the in- 
jector annuli for isothermal flow with low Mach number: 


w 


if 



( 52 ) 


The steady- state fuel injection velocity v^, which is required for the calculation of drop 
size, is computed from 


_ w, 

v f = (R g T) ,ir A f 

ci 


(53) 


Linearization and transformation of equation (52) yields 


w if (S) 




(54) 


For constant fuel flow into the fuel injector manifold, perturbations in the injector 
manifold pressure can be expressed as 


P U (S) = 



(55) 


where 

Tj steady-state turbine discharge temperature 
Vj fuel injector manifold volume 

Combining equations (54) and (55) and letting S = j<u' yield 
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, , _ _ P if (S) - P ci (S) 

Z sf (jco') = dt t + j> f = 11 w 

w if (S) 


P 2 - P 2 
Pci^f 


f P cll v /( 


j 


( 56 ) 


The imaginary part is used in the solution of equations (28) and (29). The real 
part can be compared with the critical value 9?^ resulting from the solution of 
equation (29), by neglecting the effects of changing 9t^ (change in fuel annulus area) on 
the imaginary part. 

Critical values of the pressure upstream of the injection elements and the corre- 
sponding ratios of pressure drop to chamber pressure are found from 



a ’^o 

2 


+ P 


ci 






— — 2 \l/2 

Jw, P„ 5 + Pi- 


f w rci 


ci) 



(57) 


(58) 


DIGITAL COMPUTER PROGRAM 

A digital computer program, capable of generating stability limits over a wide range 
of operating conditions and engine configurations, was developed. The stability program 
was written in FORTRAN IV language for use with the IBM 7094 computer. The program 
is structured to carry out sequentially the following operations: (1) Calculation, for a 
specified chamber pressure and from given input data, of the steady-state frequency- 
insensitive parameters required to solve equations (28) and (29), (2) calculation, for 
each selected frequency, of the feed system impedances using the procedures outlined 
in appendix C, (3) for each selected frequency, solution of equations (28) and (29) for the 
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critical parameters and 3'^, (4) translation of the parameters 3t ' Q and 9t ^ into 
pressure drop ratios using equations (57) and (58) (If the results are non-negative, com- 
pare with the steady- state operating values to gage the stability margin.) (5) repetition 
of the preceding steps for each specified chamber pressure or changes in system param- 
eters. 


Steady-State Operating Point 

Because of the unavailability of a suitable model for calculating « . , the overall 
combustion efficiency n c * or flow rates at full thrust must be given. Other required 
input data are (1) chamber, injector, and feed- system geometry, (2) propellant proper- 
ties, (3) combustion properties for the specified propellant combination, (4) oxidizer 
pump characteristic, (5) assumed fuel temperature, and (6) O/F control-valve area 
against chamber pressure for the system shown in figure 4. 

The program will calculate the steady-state frequency-insensitive parameters re- 
quired for the solution of equations (28) and (29) in the following order: 

(1) The propellant properties T , J ( 0 , and H Q are used to calculate the associated 
terms in equation (34). 

(2) The specified throat area, characteristic length L* = (V c + V n )/A^, contraction 
ratio, and total length l + l , are used to calculate the geometric term in equation (34). 

(3) For a selected chamber pressure P , the fuel temperature and O/F control- 
valve area are calculated using supplied characteristics. 

(4) For the selected chamber pressure and specified mixture ratio, c^ is com- 
puted from input data. 

(5) At full thrust, the overall efficiency 7? c+ is computed from chamber pressure 
and specified flows. At off-design chamber pressures, an assumed value of rj c * (or 
last computed value) is used to calculate flows. 

(6) From selected chamber pressure and mixture ratio, y , M , R , T , and v 

vy V/ g L L L 

are computed to allow calculation of chamber pressure at the injector face P . from 
equation (23). 

(7) Using input data for the injector, the fuel velocity v^ and fuel injector pressure 
P-j are calculated from equations (52) and (53), and oxidizer velocity and injection ori- 
fice pressure drop are calculated assuming a square-law characteristic. 

(8) Using input data for the oxidizer feed system, the steady-state pressure drops 
across the O/F control valve and dual-orifice elements (if required) are calculated. 

The program iterates to find the axial-element size which satisfies the vortex charac- 
teristic in equation (49) at the specified flow level. 
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(9) The program computes the oxidizer drop size from equation (37) using the 0. 448- 

swirler-correction factor, the burning length l from equation (41), the generalized 
length l from equation (33), the percent vaporized C from the supplied curve 

from reference 8, r? va from equation (32), and the oxidizer vaporization delay cr vo 
from equation (40). 

(10) At full thrust, T? mix is computed from rj c */v v ^p- At other chamber pressures 
a previously determined value of ?7 m ^ x is used to compute the overall efficiency, 

77 c * = ^mix^vap' ^ new va ^ ue does not agree with the assumed value in step (5) the 
new value is used to recompute flows and repeat steps (7) to (10) until agreement is 
reached. 

(11) Steady- state values of AP^ 0 /P c and AP^/P c for comparison with future com- 
puted values of the critical ratios as defined in equations (57) and (58). 

(12) The gas-residence time, 9 from equation (39) and the gas-phase mixing delay 

— o 

a m from tabular data based on figure 3 are calculated. 

(13) The engine gains X and F are computed from equations (14) and (15) with the 
partial derivative of c*^ with respect to O/F calculated by the computer at the selected 
chamber pressure and mixture ratio. 

(14) The oxidizer and fuel- system element resistances, inductances, and capaci- 
tances are calculated using the appropriate line geometry, manifold volumes, orifice, 
and valve pressure drops, etc. From the calculated oxidizer pump discharge pressure 
and specified inlet (tank) pressure, equation (48) can be used to solve for the required 
steady-, state pump speed. This speed value is then used to calculate the pump resistance 
from the slope of the supplied pump characteristic. 

(15) All pertinent information is written out. 


Calculation of Stability Limits 

After calculating all of the steady-state frequency- independent parameters re- 
quired for the solution of equations (28) and (29), a range of possible chugging frequen- 
cies to be studied must be specified. Experience has shown that a lower frequency band 
extending from 3. 5/j hertz (where r is the sum of the oxidizer -vaporization delay, 
gas-phase mixing delay, and gas residence time in milliseconds) to about 4. 0/r hertz 
and a higher band extending from 8. 7 to 12. 0/r hertz will generally enclose the solutions 
of interest. Depending on the configuration and operating point being studied, solutions 
in either or both of these banv may exist. For the studies conducted herein, the solu- 
tions within a band occurred at frequencies quite close together. For this reason a fre- 
quency increment of 0. 1 hertz was chosen in both frequency bands. 
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For the selected steady-state operating point and for- each frequency of interest, the 
following calculations were made by the computer program: 

(1) Calculate the steady- state value of the imaginary part of the fuel impedance J ^ 
using equation (56). 

(2) Calculate the real and imaginary parts of the oxidizer impedance and J 
using the reduction procedure outlined in appendix C. 

(3) Calculate the remaining frequency- dependent parameters which appear in equa- 
tions (28) and (29). 

(4) Solve equation (28) for real solutions of 9t' Q and convert to critical values of in- 
jection orifice resistance a’ using equation (51). Disregard negative values of a' 

(5) Using results from step (4), calculate the critical value (or values) of using 
equation (29). Disregard negative values of 

(6) Convert results from steps (4) and (5) into critical values of injection element 
pressure drop using equation (57) and (58). 

(7) Store the results from step (6) together with the frequency and repeat steps (1) 
to (7) for each frequency in the selected bands of frequency. 

(8) Write out the results of steps (1) to (7). An ordering subroutine is included in 
the program to write out the critical pressure-drop ratios (AP^ o /P c )' and (AP.^/P^ 5 
and the corresponding frequency in order of increasing (AP^/P )\ 

Both the steady-state and stability portions of the program are repeated for each 
chamber pressure of interest and system parameter value. The existing program is set 
up to study chamber pressures from 100 to 10 percent of the full thrust value and fuel 
annulus areas from 25 to 200 percent of the design value. 

For both the variable P c case (throttling) and the variable fuel area case, the only 
results from the stability portion of the program that are of interest are those sets of 
points ((AP ij ,/P c )', (AP io /P c )', and f = w'/27r) where (AP if /P c ) f is equal to the steady- 
state operating value corresponding to the specified fuel annulus area. This is based on 
the fact that the combustion delay times a yo and a m are sensitive to the fuel injection 
velocity and the stability limits are generated using the steady- state operating value of 
fuel velocity. For this reason, interpolation of the results written out in step (8) of the 
stability portion of the program is required. For the variable fuel area case, this inter- 
polation is built into the program. For the throttling case, the interpolation can be 
handled conveniently by hand from the output data because of the small variation of fre- 
quency along the boundary. 

Appendix D contains a list of the computer program symbols and their definitions. 
Appendix E contains a computer listing showing the actual implementation of the steps 
outlined in the steady-state and stability subprogram discussion. The indexing shown 
corresponds to the case of a dual-orifice injector configuration (see fig. 5(b)) at full 
thrust with varying fuel annulus area. Comment cards are inserted at appropriate loca- 
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tions in the listing to aid in the understanding of the program. Also shown is a printout 
of the resultant stability limits for this system together with other pertinent informa- 
tion. The computer program as shown in the listing, is set up for the expander cycle 
with gaseous fuel being fed into the injector manifold at a constant rate. The program- 
ming of the general liquid-liquid system, as described in the analysis, requires the 
calculation of both oxidizer and fuel vaporization delays, gas-phase delays, feed system 
impedances, etc. 

Appendix F contains a flow chart for the digital program and a discussion, which is 
intended as a guide for users of the pregram. This information, together with the com- 
ment cards contained in the listing (appendix E), should help the user to modify the pro- 
gram to study different engine configurations. The steps required to analyze liquid- 
liquid propellant systems are also discussed in appendix F. 


RESULTS AND DISCUSSION 
Single-Orifice Configuration 

A flox-methane engine system using an expander cycle and an oxidizer injector as 
shown in figure 5(a) was studied using the digital program. Table I shows the specified 
full thrust values of the chamber pressure, flows, etc. 

At full thrust (P c = 3. 45x10® N/m^) the steady- state and stability portions of the 
program were used to determine critical injector pressure drops for a range of fuel in- 
jector annulus areas. Figure 6 is a plot of the digital data. For the nominal fuel area, 
the ratio of fuel- injector-pressure drop to chamber pressure (AP.j/P c ) is 0. 152 at full 
thrust with a nominal oxidizer ratio of 0. 250. For the nominal fuel area, the stability 


TABLE I. - SINGLE-ORIFICE CONFIGURATION OPERATING CONDITIONS AT FULL THRUST 


Q Z* 

Chamber pressure, P c , N/m 3. 45x10 

Mixture ratio, O/F 5. 25 

Oxidizer mass flow rate, w Q , kg/sec 4. 68 

Fuel mass flow rate, w f , kg/sec 0.891 

I o c 

Fuel-injector total pressure, P^, N/m 4. 02x10 

Fuel-injector temperature, T, f , K 763 

II n z? 

Oxidizer-injector total pressure, P. , N/m^ 4.35x10 

lO n a 

Oxidizer pump discharge pressure, P Q p d , N/m 4. 69x10 

Oxidizer tank pressure, P To , N/m 2 2. 76x10^ 

Oxidizer pump speed, N Q , rpm 2. 465x10^ 

Oxidizer velocity, v , m/sec 20. 42 

Fuel velocity, Vp m/sec 214.4 
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Figure 6. - Comparison of analog and digital stability limits for single- 
orifice, full-thrust injector with varying fuel area. Delays sensitive 
to fuel area. 


program yielded a critical oxidizer ratio of 0. 130, which indicates stable operation with 
the nominal geometry. At the critical oxidizer ratio, a chugging frequency of 210. 5 
hertz was determined by the program. Changes in the oxidizer pressure drop are as- 
sumed to be made by changing the size of the oxidizer injection orifice without affecting 
the injection velocities, delay times, etc. If one assumes that the resulting changes in 
pressures, upstream of the injector, have only second-order effects on the stability of 
the engine, then the difference between the operating point and critical pressure drops 
represents the stability margin at the operating point. A stability boundary was formed 
by computing critical oxidizer pressure drops for various fuel annulus areas. For each 
selected fuel area, the delay times were calculated using the steady-state operating 
point portion of the program. Figure 6 also shows a higher frequency boundary at lower 
oxidizer and fuel pressure-drop ratios. For these operating conditions, the higher fre- 
quency boundary is not important since the lower frequency boundary will determine the 
stability of operation. For comparison with digital results, results from an analog com- 
puter simulation are also shown. On the analog computer, a small-impulse disturbance 
was introduced in the rate of change of burned products (Aw^ as shown in fig. 2) with a 
resulting sinusoidal oscillation in chamber pressure. The oxidizer injection orifice 
size was slowly decreased while adjusting the pump speed to maintain flow until steady- 
amplitude oscillations in chamber pressure were observed. Frequencies on the analog 
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were determined by using an electronic counter which gave a visual indication of the 
oscillation period (in msec). Excellent agreement in the boundary location and frequency 
is indicated in figure 6, which substantiates the inference of stability margin from the 
digital data. In terms of stability, a slight advantage in decreasing the fuel-injector 
pressure drop is indicated. A study of the computer results, however, indicated a loss 
in combustion efficiency due to a decrease in the injection momentum ratio. This re- 
sults in an increase in drop-size, and in decreases in the generalized length parameter 
l and fraction vaporized 0 for the oxidizer. Although the combustion efficiency 

& t;Il 

due to mixing is assumed constant (r/ m ^ x = 0. 973), the decrease in momentum ratio 
would probably result in a further decrease in performance due to incomplete mixing. 

The effect of throttling the engine at constant mixture ratio was also studied. A 
turbine-bypass valve (see fig. 4) was assumed to control chamber pressure through 
pump speed. A mixture ratio control valve (fig. 5(a)) area was varied in a specified 
manner to maintain the chamber mixture ratio. Using the nominal fuel area, the steady- 
state and stability portions of the program were run for various chamber pressures. 
Figure 7 shows the results from the digital program. Figure 7(a) shows a plot of the 
critical values (AP^/P c )' and (AP io /P c )' full thrust together with the frequency 
f' = w'/2v. Since all points on the boundary were computed using the operating point 
delay times, the only boundary point of interest corresponds to the operating value of 
fuel-injector pressure drop. 

Interpolation on the plot yields the critical value of oxidizer pressure-drop ratio 
equal to 0. 130 and frequency about 210 hertz which agrees with the results shown in 
figure 6. As chamber pressure was decreased (figs. 7(b) and (c)) a shift in the lower 
frequency boundary to the left is observed with the operating point moving upward and to 
the left. The upward movement of the operating point is due to the assumed static rela- 
tion between turbine discharge temperature and chamber pressure. At a throttling ratio 
of 3. 33, figure 7(d) shows the appearance of a higher frequency boundary from the left. 

No lower frequency boundary is indicated since at this operating point, the valve and 
pump impedances are sufficient to stabilize the lower frequency mode. Figure 7(e) 
shows that at a throttling ratio of 5, a higher frequency instability is predicted at a fre- 
quency of 529 hertz. Figure 7(f) is a plot of the critical oxidizer pressure drop ratios 
at the operating point value of fuel pressure drop against the throttling ratio. By ne- 
glecting the second-order effects on stability of changing the pressure level upstream of 
the injector, the difference between the throttling path and the boundaries can be con- 
sidered to be the stability margin of the engine. For throttling ratios between 2.0 and 2. 5, 
no injection orifice pressure drop is required. The limiting factor on throttling appears 
to be the higher frequency instability with a throttling limit of 4. 82. 


27 



Critical ratio of oxidizer injection orifice pressure to chamber pressure Ratio of fuel -element pressure change to chamber pressure 
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(a) Throttling ratio, 3.45xl0 6 /P c = 1.0. (b) Throttling ratio, 3.45xl0 6 /P c = 1.25. 


*c) (c) Thruster ratio, 3.45xlofyp c - 1.67. 
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(d) Throttling ratio, 3.45xKr/P c = 3.33. 


(e) Throttling ratio, 3.45xlO b /P c = 5.0. 



(f) Throttling limits for single-orifice configuration (throttling ratio = 3.45xl0 6 /P c ). 


Figure 7. - Digital results for stability limits for single -orifice, nominal -fuel area injector. Delays evaluated at operating point. 
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Dual -Orifice Configuration 


A flox-methane engine system using an expander cycle and an oxidizer injector as 
shown in figure 5(b) was studied using the digital program. Table II shows the specified 
full thrust operating conditions for this configuration. Figure 8 shows the plot of digital 
data showing the effect of varying fuel annulus area and oxidizer injection orifice size on 
the stability limits at full thrust. For this configuration more advantage is indicated for 
decreasing the fuel injector pressure drop ratio. Digital data indicated that operation at 
a fuel injector pressure drop ratio of 0. 15 (same as single- orifice case) could be accom- 
plished with no performance degradation due to incomplete vaporization. Appendix E 
contains a printout of the digital results for this case. Figure 8 also shows the higher 
frequency boundary dominating the stability at fuel-injector pressure drops below about 
0. 1 P . Analog computer results are also shown in figure 8. Excellent agreement be- 
tween analog and digital data were obtained. Unfortunately, a limitation of the analog time 
delay capability prevented any determination of the higher frequency boundary on the ana- 
log. Analog computer operation in the stable-low, unstable- high frequency region has 
been achieved for other engine configurations. Figure 9 shows the response of chamber 
pressure to an impulse distance in flow for oxidizer injector pressure-drop ratios both 
greater than and less than the critical value. Figure 9(a) shows the decay of oscillations 
for stable operation with both frequencies apparent in the trace. Figure 9(b) shows a 
decay in the low-frequency component but a growing higher frequency component. 

The effect of throttling the dual orifice system at constant mixture ratio was also 
studied using the digital program. As in the case of the single orifice injector a turbine 
bypass throttle was assumed with mixture-ratio controlled by a valve in the secondary 


TABLE n. - DUAL ORIFICE CONFIGURATION OPERATING CONDITIONS AT FULL THRUST 


9 

Chamber pressure, P c , N/m 

Mixture ratio, O/F 

Oxidizer mass flow rate, w , kg/sec 

Fuel mass flow rate, w^, kg/sec 

Primary flow mass flow rate, w , kg/sec . . . 

Secondary flow mass flow rate, w 0 . kg/sec . . 

sec n 

Secondary total pressure, P , N/rn 

s ec n 

Oxidizer injector total pressure, P. , N/m^ . . 

lO c\ 

Fuel injector total pressure, P.^, N/irr 

Fuel injector temperature, T.,, K 

Oxidizer pump discharge pressure, P ( j, N/m 

Oxidizer tank pressure, P To , N/m 2 

Oxidizer pump speed, N Q , rpm 

Oxidizer velocity, v Q , m/sec 

Fuel velocity, v^, m/sec 


3.46X10 6 
. . 5.25 

. . 4. 74 

. . 0.903 
. . 0.479 
. . 4.27 

. . 4.03 

3.84X10 6 
4. 52X10 6 
... 747 

4. 35X10 6 
2. 76X10 5 
2. 396X10 4 
. . 12.92 
. . 217. 1 
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Ratio of fuel element pressure change to chamber pressure 


Digital data at critical frequency (f 1 ) 
O Analog data at analog frequency (f g ) 
□ Nominal operating point 


Frequencies, 



Ratio of oxidizer injection orifice pressure change 
to chamber pressure 


Figure 8. - Comparison of analog and digital stability 
limits for dual-orifice, full-thrust injector with vary- 
ing fuel area. Delays sensitive to fuel area. 



(a) Oxidizer injection orifice pressure change ireaterthan 
critical. 



(b) Oxidizer injection orifice pressure change less than 
Critical. 

Figure 9. - Response of chamber pressure to pulse-type 
disturbance on analog for conditions where higher fre- 
quency is unstable and lower frequency is stable. 


(axial) flow path (see fig. 5(b)). Figure 10 shows the results of the throttling study. 
Figures 10(a) and (b) show stable operation for throttling ratios of 1 to 1. 43. At a 
throttling ratio of 3. 33, the higher frequency boundary appears from the left in figure 
10(c) with frequencies around 502 hertz at the boundary. Figure 10(d) shows near 
neutrally stable operation at a throttling ratio of 5. 0 with a chugging frequency of 
493 hertz. Figure 10(e) shows the movement of the operating point and stability margin 
for varying thrust. Figure 10(e) indicates a throttling limit of 4. 8 which agrees with the 
single- orifice result. 
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(b) Throttling ratio, 3.45xl0 6 / (c) Throttling ratio, 3.45xl0 6 / (d) Throttling ratio, 3. 45xl0 6 / 

P c - 1. 43. P c = 3. 33. P c “ 5.0. 



Throttling ratio 

(e) Throttling limits for dual-orifice configuration (throttling ratio s 3.45xl0 6 /P c ). 

Figure 10. - Digital results for stability limits for dual-orifice, nominal-fuel area injector. Delays evaluated at operating point. 
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CONCLUDING REMARKS 


A digital program has been developed to determine stability limits for throttleable, 
bipropellant rocket engines. The program was written in FORTRAN IV language for 
use with the IBM 7094 computer. 

The program consists of (1) a steady- state operating point subprogram which com- 
putes required parameters for solution of the characteristic equation and (2) a stability 
subprogram which solves the characteristic equation for critical values of the injection 
element pressure drops over a range of chugging frequency. Guidelines are provided 
for choosing the range of frequency to be studied. 

Available vaporization and drop size correlations are used to compute vaporization 
efficiencies and combustion delay times. The effects of injector geometry, valve re- 
sistance, pump characteristics, suction line inductance, etc. , are included in the basic 
program. 

Digital stability data were computed for an expander cycle engine using (1) a con- 
ventional single-orifice concentric -tube injector and (2) a dual-orifice concentric-tube 
injector. The effects of throttling at constant mixture ratio and of varying the injector 
pressure drops at full thrust on stability were determined. Excellent agreement between 
digital and analog computer data were obtained. 

For both injector configurations, throttling was limited to about 5 because of higher 
frequency instabilities (about 500 Hz), which is characteristic of the double -dead- time 
model used in the program. For the dual- orifice configuration, results indicated that 
the design-point fuel-injector pressure drop could be decreased by 50 percent without 
degrading performance or stability. 

The development and workings of the computer program were discussed in detail 
and a program listing, together with a typical printout of results are provided. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, September 8, 1970, 

128-31. 
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APPENDIX A 


o 

A cross-sectional area, m 

real part of oxidizer feed system 
impedance, excluding injection 
orifice resistance, 
((N)(sec)/(m 2 )(kg)) 

2 

Aj. throat area, m 

a oxidizer injection orifice resis- 

tance, ((N)(sec)/(m 2 )(kg)) 

B bulk modulus of liquid, N/m 2 

& imaginary part of oxidizer feed 

system impedance, excluding 
injection tube reactance, 
((N)(sec)/(m 2 )(kg)) 

b oxidizer injection tube inductance, 

((N)(sec 2 )/(m 2 )(kg)) 

c* characteristic velocity, m/sec 

d diameter, m 

F partial derivative of chamber pres- 

sure with respect to fuel flow as 
defined by eq. (15), 
((N)(sec)/(m 2 )(kg)) 

fraction of fuel vaporized 

f frequency, Hz 

denotes functional relation 

g gravitational acceleration, 

9. 8 m/sec 2 

g c gravitational constant, 

1 (kg)(m)/(N)(sec 2 ) 

H heat of vaporization for liquid, J/kg 


imaginary part of feed system im- 
pedance, ((N)(sec)/(m 2 )(kg)) 

chamber momentum loss coeffi- 
cient as defined by eq. (23) 

Kj fuel-injector flow coefficient as 

defined by eq. (52), 
((kg)(m 2 )(K 1 / 2 )/(sec)(N)) 

K r drop size coefficient as defined by 
eq. (35) 

K^p pump pressure rise coefficient as 

defined by eq. (48), 
((N)(sec 2 )/(m 2 )) 

Kgp pump flow coefficient as defined by 
eq. (48), kg' 1 

L* characteristic length, m 

l length, m 

j( molecular weight for liquid, 

((kg)/(kg)(mole)) 

M Mach number 

N pump speed, sec" 1 

C fraction of oxidizer vaporized 

O/F mixture ratio 

2 

P total pressure, N/m 

p static pressure, N/m 2 

Rg gas constant, J/(kg)(K) 

r drop radius, m 

dt real part of feed system impe- 

dance, ((N)(sec)/(m 2 )(kg)) 

S Laplace operator, sec" 1 


SYMBOLS 

J 
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Sf nozzle shape factor = V n /A c £ n 

T temperature, K 

t time, sec 

3 

V volume, m 

v velocity, m/sec 

W mass, kg 

w mass flow rate, kg/sec 

X partial derivative of chamber pres- 

sure with respect to oxidizer flow 
as defined by eq. (14), 
((N)(sec)/(m 2 )(kg)) 

Y substitution parameter as defined 

in eq. (5) 

Z complex flow impedance, 

((N)(sec)/(m 2 )(kg)) 

a real part of general system element 

o 

impedance, ((N)(sec)/(m )(kg)) 

/3 imaginary part of general system 

element impedance, 
((N)(sec)/(m 2 )(kg)) 

y specific heat ratio 

AP total pressure loss (rise for pump), 

N/m 2 

Aw disturbance in mass flow rate, 
kg/sec 

6 surface tension of liquid, N/m 

e c contraction ratio 

t] efficiency 

gas residence time as defined by 
eq. (34), sec 

X real part of complex variable S, 
rad/sec 

2 

ju viscosity of liquid, ((N)(sec)/(m )) 


p 

o 

mass density, kg/m 

a 

time delay, sec 

T 

sum of oxidizer delays and gas 
residence time, sec 

<P 

flow parameter as defined by 
eq. (43) 

* 

pressure- rise parameter as de 
fined by eq. (43) 

u> 

angular frequency, rad/sec 

Subscripts: 

a 

analog value 

b 

burned 

c 

chamber (at nozzle inlet) 

ci 

chamber (at injector face) 

cr 

critical 

c* 

characteristic velocity 

f 

fuel 

fb 

burned fuel 

gen 

generalized 


indices for general feed system 
elements 

if 

fuel injector 

io 

oxidizer injector 

ip 

injected propellant droplet 

m 

gas-phase mixing and reaction 

man 

manifold 

mix 

due to incomplete mixing 

n 

nozzle 

0 

oxidizer 

°2 

liquid oxygen 

ob 

burned oxidizer 
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opd 

oxidizer pump discharge 

V 

vaporized 

P 

pump 

vap 

due to incomplete vaporization 

pr 

primary flow 

vf 

vaporized fuel 

s 

suction line 

VO 

vaporized oxidizer 

sec 

secondary flow 

50 

50 -percent vaporized 

sf 

fuel system output 

Superscripts: 

so 

oxidizer system output 

- 

mean or steady- state value of any 

T 

total 


variable 

th 

theoretical 

T 

critical value of any variable 

V 

valve 


time-varying portion of any vari- 


able (i. e. , x = x + x) 
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APPENDIX B 


SOLUTION OF CHARACTERISTIC EQUATION 


The characteristic equation describing the bipropellant engine system was derived 
in the text and was 


-1 = 


XK c -< c 'vo«m> S FK c 
e + e 


z so ® 


Z sf (s) 


e e s + 1 


(Bl) 


By letting S = JW, <Vo = "'vo + %• a Tf “ a vf + “V Z so = «o ' iJ o< and Z sf=*t + H- 
the characteristic equation becomes 




XK C cos(wa To ) - jXK c sin(uxr To ) FK C cos(wa Tf ) - jFK c sin(a)a Tf ) 


*o + j y o 




(B2) 


By multiplying both sides of equation (B2) by the terms (&' + and sepa- 

rating both sides into real and imaginary parts, and equating reals and imaginaries, the 
real part of the fuel impedance, 9t'^ can be solved for using the real and imaginary equa- 
tions. From the imaginary equation, 


+ ^o ' XK c sin ( a) ' a To^ 


St\ FK c ^’ 0 sin(co’a Tf ) - FK^ cos(w’CT Tf ) - XK^ cos(u)'a To ) + a>»e S Q J t - 


(B3) 


From the real equation, 


co 


Vo - *o - XK c 


FK c ^' q cos(w’u Tf ) + FK c jr Q sin(w’ff Tf ) + XK c > f sin(w'cr To ) -S Q S f - <*>'8 g S f & 0 

(B4) 
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Equations (B3) and (B4) yield the following quadratic equation in terms of 3' q : 



where 


j 

FK, 


l (l + o> ,2 0 2 )j + ^|-XK C sin(a)'a To - u>'a T{ ) - _^J[cos(w'a To ) - w’0 g sin(u> , CT To jjj. 

{ _ 2x7 > { r _ _ -i 

- > 0 XK C cos(a)’a To - co'<J Tf ) + =2 — |sin(a)'cr To ) + w'0 g cos(a>’ff To )J 


J y^ 

Yo 

fit 



X 2 ^,K 


f C l 

- J 


(B5) 


= o 


Y = c o'Q cos(co’CTrpj) + sinCco'CTrpj) 

D 
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APPENDIX C 


EVALUATION OF FEED SYSTEM IMPEDANCES 

Solution of the characteristic equation for the bipropellant system, as outlined in 
appendix B, requires the knowledge of the imaginary parts of the feed system imped- 
ances and at each frequency of interest. The following discussion will outline 
the steps required to evaluate these parameters for both the single and dual injector 
configurations. In general, the evaxuation will involve (1) the breaking up of the feed 
system in question into its elements, each having a flow impedance Z. with the real 
and imaginary parts a ^ and /3p respectively, (2) the manipulation of these elements 
(if necessary) into series and parallel combinations, and (3) the stepwise reduction of 
these combinations to obtain the feed-system impedance (both real and imaginary parts) 
at the frequency of interest. 


Injector 

elements 


0/F control 
valve 


Pump 


Suction 

line 

z = a + jo;b 


_ 

z v 

_ 

_ 

z p 


2 S 


Injector 


Pump 

manifold 


discharge 

volume 


volume 

z man 


z opd 


Figure 11. - Impedance representation of single orifice oxidizer feed system. 


The step-wise reduction of the system shown in figure 11 would be 


Z 1 = Z p + Z s 

(Cl) 

z _ Z l Z opd 

(C2) 

2 Z l +Z opd 


Z 3 = Z 2 + Z V 

(C3) 
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<7 Z 3 Z man 

z, 4 

Z- + Z mQ¥1 
o man 

Z so = Z 4 + a + j wb (C 5) 

For this case, the imaginary and real parts of Z would be 

SO 

J Q = cob + 0 4 = a + a 4 (C6) 

where /3 4 is the imaginary part of the impedance Z^, evaluated at the frequency of in- 
terest. To evaluate the real and imaginary parts of the impedances, the following equa- 
tions are very useful: If 


S' 

(C4) 


z k = z i + Z j 


then 




Q-i + 


0/i 




If 


Then 


. (»i ♦ 0?) + »i (»? + »j) 

O'k “ ' 

(<*1 + a-j) 2 + Oj + /3j) 2 

+ (Vj) 2 + (/3j + /3j) 2 


(C7) 


(C8) 


Equations (C7) and (C8) form the basis for the feed- system impedance calculations at the 
frequency of interest regardless of the feed system configuration. 
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Figure 12. - Impedance representation of dual-orifice oxidizer feed system. 



Figure 13. - Modified impedance representation of dual -orifice oxidizer feed system. 
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Figure 12 shows the elements for the dual-orifice oxidizer feed system. Since the 
stepwise reduction of the feed- system elements is based on parallel and series combina- 
tions, the tee, formed by the secondary elements, manifold, and control valve must be 
converted to a delta form. Figure 13 shows the modified feed system representation for 
the dual orifice system. 

The impedances, Zj, Zg, and Z 3 , are evaluated from the following conversion 
equations: 


Z i — z + z + 
1 sec man 


Z Z 
sec man 


J V 


(C9) 


Z 2 - Z v + Z sec + 


Z Ztt 
sec V 


man 


(CIO) 


Zg - Zy + Z man + 


Z V Z man 


sec 


(Cll) 


The stepwise reduction of the modified system is 


Z 4 - V Z S 


Z 5 = 


Z 4 Z opd 
Z 4 + Z opd 


Z 6 = 


Z 5 Z 3 
Z 5 + Z 3 




Z 2 + Z pr 


Z g - Z g + z 7 


Z 9 = 


Z 1 Z 8 
Z 1 + Z 8 


(C12) 

(C13) 

(C14) 

(C15) 

(C16) 

(C17) 
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Z so - z 9 + a + 3 wb 


(C 18 ) 


y Q = cub + 0 9 (C19) 

~ a + 

As discussed with the single- orifice configuration, equations (C7) and (C8) were used to 
evaluate the various a^ T s and /^'s at each selected frequency. 

After solutions to the characteristic equation have been found, the critical values of 
the real part of impedance &' Q are converted to the corresponding critical values of the 
injection-element resistance a’ using either equation (C6) or (C19). 
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APPENDIX D 


AC 

ACCH 

ACCL 

ACOUS 

ACRH 

ACRL 

AG 

AHX 

AKX 

ALPHA 

ALPHAG 

AOV 

AOV1 

APF 

APFO 

AT 

AX 

AXA 

AYX 

A1 

A2 


COMPUTER SYMBOLS 

2 2 

combustion chamber cross-sectional area, in. (m ) 

storage location for critical oxidizer pressure drop ratio at operating fuel 
pressure-drop ratio (high frequency range) 

storage location for critical oxidizer pressure drop ratio at operating fuel 
pressure -drop ratio (low frequency range) 

acoustic velocity at nozzle entrance, in. /sec (m/sec) 

storage location for critical oxidizer pressure drop ratio (high frequency 
range) 

storage location for critical oxidizer pressure-drop ratio (low frequency 
range) 

coefficient of second- order term in eq. (28) 

o 

imaginary part of dual orifice Zg in appendix C, sec/in. 
((N)(sec)/(m 2 )(kg)) 

o 

imaginary part of dual orifice Z^ in appendix C, sec/in. 
((N)(sec)/(m 2 )(kg)) 

factor under square root in eq. (28) 

product of angular frequency and gas residence time, rad 

2 2 

control valve area, in. (m ) 

2 2 

full thrust control-valve area, in. (m ) 

2 2 

fuel injector annulus area per element, in. (m ) 

2 2 

design value of fuel annulus area per element, in. (m ) 

2 2 

nozzle throat area, in. (m ) 

real part of single- orifice Z^ or real part of dual- orifice Z^ in appen- 
dix C, sec/in. 2 ((N)(sec)/(m 2 )(kg)) 

2 2 

real part of oxidizer impedance, sec/in. ((N)(sec)/(m )(kg)) 

real part of dual-orifice Z „ in appendix C, sec/in. ((N)(sec)/(m )(kg)) 

cosine of phase angle due to gas-phase mixing delay 

cosine of phase angle due to oxidizer vaporization delay 
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A3 

BG 

BO 

BX 

BXA 

BXI 

B1 

B2 

B3 

CDF 

CDO 

CDP 

CDS 

CG 

COEFF 

COF 

COM 

COP 

CSTAA 

CSTAR 

CX 

Cl 

C2 

DDX 

DENE 

DO 


cosine of phase angle due to total oxidizer delay 
coefficient of first order term in eq. (28) 
oxidizer bulk modulus, (N/m ) 

2 

imaginary part of single -or if ice Z^ in appendix C, sec/in. 
((N)(sec)/(m 2 )(kg)) 

2 2 

imaginary part of oxidizer impedance, sec/in. ((N)(sec)/(m )(kg)) 

2 

reciprocal of imaginary part of dual orifice Z^ in appendix C, in. /sec 
( (m 2 ) (kg) /(N) ( s ec) ) 

sine of phase angle due to gas -phase mixing delay 
sine of phase angle due to oxidizer vaporization delay 
sine of phase angle due to total oxidizer delay 
flow coefficient of fuel-injector elements 
flow coefficient of oxidizer injection orifices 
flow coefficient of primary injector elements 
flow coefficient of secondary injector elements 
coefficient of zero- order term in eq. (28) 
gas-phase delay coefficient, sec 
fuel-injector manifold capacitance, in. (m ) 

single- orifice injector manifold or dual-orifice secondary manifold capaci- 
2 2 

tance, in. (m ) 

single -orifice oxidizer pump discharge or dual-orifice pump discharge plus 

2 2 

primary manifold capacitance, in. (m ) 
storage location for characteristic velocity, ft/sec (m/sec) 
characteristic velocity, ft/sec (m/sec) 

O 

real part of single- orifice Zg or dual- orifice Zg in appendix C, sec/in. 
((N)(sec)/(m 2 )(kg)) 

oxidizer pump flow coefficient 

oxidizer pump head rise coefficient 

denominator of dual orifice impedance Zg in appendix C 
denominator of efficiency due to incomplete vaporization 
oxidizer injection orifice diameter, in. (m) 
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DP 

DPF 

DPFC 

DPF1 

DPF2 

DPJ 

DPO 

DPOC2, 5, 8 

DPP 

DPS 

DS 

DSU 

DX 

EC 

ED 

EFF 

EFF1 

ELC 

ELCT 

ELEFF 

ELGEE 

ELJ 

ELLS 

ELN 

ELO 

ELSU 


primary injector element diameter, in. (m) 

ratio of fuel injector pressure drop to chamber pressure at throat 

critical value of fuel-injector pressure drop to chamber pressure ratio 

critical value of fuel-injector pressure drop to chamber pressure ratio 
when there is more than one solution 

critical value of fuel-injector pressure drop to chamber pressure ratio 
when more than one solution 

9 

oxidizer injection- orifice pressure drop, psi (N/m ) 

ratio of oxidizer injection- orifice pressure drop to chamber pressure 

ratios of oxidizer injection- orifice pressure drop to chamber pressure 
when more than one solution 

2 

pressure drop from oxidizer pump discharge to injector face, psi (N/m ) 

pressure drop from oxidizer pump discharge to injection orifice input, 
psi (N/m 2) 

secondary injector- element diameter, in. (m) 
suction line diameter, in. (m) 

imaginary part of single orifice Zg or dual orifice Zg in appendix L 
sec/in. 2 ((N)(sec)/(kg)(m 2 )) 

contraction ratio 

element density 

combustion efficiency 

storage location for combustion efficiency 
cylindrical chamber length, in. (m) 

total chamber length from injector face to throat, in. (m) 
effective length (ref. 7), in. (m) 
generalized length (ref. 7), in. (m) 

2 2 2 2 

oxidizer injection- element inductance, sec /in. (sec /m ) 

suction line length, in. (m) 

conical nozzle length, in. (m) 

oxidizer injection-element length, in. (m) 

suction line inductance, sec 2 /in. 2 (sec 2 /m 2 ) 
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ELSTR 

ELT 

ELX 

EL50 

EMAGO 

EMAGF 

EMC 

EMR 

EMXI 

ENE 

ENF 

ENO 

ENUM 

ENX 

ETAM 

ETAV 

EX 

FFH 

FFL 

FG 

FH 

FL 

FOP 


characteristic length, in. (m) 

chamber-nozzle geometry coefficient (ref. 8), in. (m) 

o 

real part of dual orifice impedance Zg in appendix C, sec/in. 

- ((N)(sec)/(m 2 )(kg)) 

*' length required to vaporize 50 percent of oxidizer droplet mass, in. (m) 

2 2 

imaginary part of oxidizer impedance, sec/in. ((N)(sec)/(m )(kg)) 

2 2 

imaginary part of fuel impedance, sec/in. ((N)(sec)/(m )(kg)) 

Mach number at nozzle entrance 
injectffr momentum ratio 

2 

reciprocal of imaginary part of dual orifice Zg (see appendix C), in. /sec 
((m 2 )(kg)/(N)(sec)) 

number of injector elements 

number of fuel elements 

number of oxidizer elements 

numerator of efficiency due to incomplete vaporization 

2 

real part of dual-orifice impedance Z^ (see appendix C), sec/in. 
((N)(sec)/(m 2 )(kg)) 

efficiency due to incomplete mixing 

efficiency due to incomplete vaporization 

2 

real part of single-orifice Z2 or dual-orifice Zg (see appendix C), sec/in. 
((N)(sec)/(m 2 )(kg)) 

storage location for frequency at critical oxidizer pressure-drop ratio, oper- 
ating fuel pressure-drop ratio (high-frequency range), Hz 

storage location for frequency at critical oxidizer pressure-drop ratio, oper- 
ating fuel pressure-drop ratio (low-frequency range), Hz 

2 

partial derivative of chamber pressure with respect to fuel flow, sec/in. 
((N)(sec)/(m 2 )(kg)) 

storage location for frequency at critical injector pressure-drop ratios (high- 
frequency range), Hz 

storage location for frequency at critical injector pressure-drop ratios (low- 
frequency range), Hz 

fraction of oxidizer flow through primary (tangential) injector elements 
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FREQ 

FX 

GAM 

GAMMA 

G 

GC 

GX 

HV 

HVT 

I 

II 
IL 
IM 
IX 
IX 1 
IZ 
J 

JJ 

JM 

JN 

JX 

K 

KK 

KX1 

L 

LI 

LJ 

LX 


frequency, Hz 

imaginary part of single orifice Z 2 or dual orifice Z g (see appendix C), 
sec/in. 2 ((N)(sec)/(m 2 )(kg)) 

input value of specific heat ratio for combustion products 
specific heat ratio of combustion products 
acceleration due to gravity, ft/sec (m/sec) 

gravitational constant, 32. 2 ((lbm)(ft)/(lbf)(sec 2 )) (1 (kg)(m)/(N)(sec 2 )) 

2 

real part of dual orifice impedance Zg (see appendix C), sec/in. 
((N)(sec)/(m 2 )(kg)) 

oxidizer heat of vaporization, Btu/lbm (J/kg) 
heat of vaporization coefficient (ref. 8) 
throttling index 

multiple oxidizer solution index 
high-frequency solution ordering index 
multiple solution frequency index 
low-frequency solution ordering index 
low-frequency solution order index 
variable fuel-area solution order index 
calculation counting index 
type of study index 

existance of high-frequency solution index 
existance of low-frequency solution index 
low-frequency solution ordering index 
variable fuel-area index 
high-frequency solution ordering index 
high-frequency solution ordering index 
solution ordering index 
low-frequency solution ordering index 
high-frequency solution ordering index 
high-frequency solution ordering index 



MM 


type of injector index 
MN lower throttling limit index 

N number of boundary points 

ND number of solution points in set 

NE one more than number of boundary points 

NN frequency band index 

NS one less than number of low-frequency boundary points 

Nl number of low-frequency boundary points 

N2 number of high-frequency boundary points 

Nil one less than number of low-frequency boundary points 

N22 one less than number of high-frequency boundary points 

OF mixture ratio 

OX imaginary part of dual orifice Zg in appendix C, sec/in. ((N) (sec) / (kg) (m ) ) 

2 

PC total chamber pressure at nozzle entrance, psi (N/m ) 

2 

PCI total chamber pressure at injector face, psi (N/m ) 

o 

PCI full thrust total chamber pressure at nozzle entrance, psi (N/m ) 

2 

PF fuel-injector manifold pressure, psi (N/m ) 

2 

PO oxidizer pump discharge pressure, psi (N/m ) 

2 

PS oxidizer pump inlet pressure, psi (N/m ) 

PSUP pressure feeding secondary (ax'-al) injector elements, psi (N/m ) 

PV percent vaporized with chamber, nozzle 

2 2 

PX real part of dual orifice Zg (see appendix C), sec/ in. ((N)(sec)/(kg)(m )) 

^ 2 

QX imaginary part of dual orifice Zg (see appendix C), sec/in. 

((N)(sec)/(kg)(m 2 )) 

RC combustion products gas constant, ft/°R (m/K) 

RESF linearized fuel 

2 2 

RESOB linearized secondary-element resistance, sec/in. ((N)(sec)/(kg)(m )) 

2 2 

RESOJ linearized oxidizer injection-orifice resistance, sec/in. ((N)(sec)/(kg)(m )) 

2 2 

RESOP linearized oxidizer pump resistance, sec/in. ((N)(sec)/(kg)(m )) 
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RE SOS 

RESOV 

RF 

RHO 

RMCTN 

RPC 

RPC 1,2 

RRF 

RRF1, 2 

RX 

RX1, 2 
S 

SLOPA 

SLOPE 

SPEED 

ST 

SWC 

TAUM 

TAUT 

TAUV 

TC 

TEMP 

TF 


linearized single- orifice valve resistance or dual- orifice primary element 

2 2 

resistance, sec/in. ((N)(sec)/(kg)(m )) 

2 2 

linearized dual-orifice valve resistance, sec/in. ((N)(sec)/(kg)(m )) 
fuel gas constant, ft/°R (m/K) 

O O 

oxidizer mass density, lbm/in. (kg/m ) 
mean oxidizer drop radius, in. (m) 

critical value of real part of oxidizer injection-orifice impedance (one so- 
lution), sec/in. 2 ((N)(sec)/(kg)(m 2 )) 

critical values of real part of oxidizer injection- orifice impedance (mul- 

2 2 

ti pie solutions), sec/in. ((N)(sec)/(kg)(m )) 

critical value of real part of fuel-injector impedance (one solution), 
sec/in. 2 ((N)(sec)/(kg)(m 2 )) 

critical values real part of fuel-injector impedance (multiple solutions), 
sec/in. 2 ((N)(sec)/(kg)(m 2 )) 

critical value of real part of oxidizer feed- system impedance (one solution), 
sec/in. 2 ((N)(sec)/(kg)(m 2 )) 

critical values of real part of oxidizer feed- system impedance (multiple 
solutions), sec/in. 2 ((N)(sec)/(kg)(m 2 )) 

nozzle shape factor (ref. 8) 

storage location for partial derivative of characteristic velocity with re- 
spect to mixture ratio, ft/sec (m/sec) 

partial derivative of characteristic velocity with respect to mixture ratio, 
ft/sec (m/sec) 

oxidizer pump speed, rpm 

oxidizer surface tension, dyne/cm (N/m) 

correction factor for swirler 

gas-phase mixing and reaction time delay, sec 

total oxidizer time delay, sec 

oxidizer vaporization time delay, sec 

combustion temperature, °R (K) 

temporary storage location for ordering solutions 

fuel temperature, °R (K) 
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TFl 

THETAG 

THETAM 

THETAV 

TK 

TPE 

TT 

UF 

UPH 

UPL 

UPPH 

UPPL 

VC 

VELC 

VF 

VIS 

VN 

VOF 

VOM 

VORTC 

VOS 

w 

WF 

WFl 

WO 


full thrust fuel temperature, °R (K) 

gas residence time, sec 

angle due to gas -phase mixing delay, rad 

angle due to oxidizer vaporization delay, rad 

input value of combustion temperature, K 

thrust per element, Ibf (N) 

temperature coefficient (ref. 8) 

unit correction factor, 12 in. /ft (1 m/m) 

storage location for high-frequency solution of critical fuel-injector pres- 
sure drop- ratio 

storage location for low-frequency solution of critical fuel-injector 
pressure-drop ratio 

storage location for ordered values of operating fuel-injector pressure- 
drop ratio (high-frequency range) 

storage location for ordered values of operating fuel-injector pressure- 
drop ratio (low-frequency range) 

3 3 

cylindrical chamber volume, in. (m ) 

gas velocity at nozzle entrance, in. /sec (m/sec) 

fuel-injection velocity, in. /sec (m/sec) 

oxidizer viscosity, cp ((N)(sec)/m ) 

3 3 

conical nozzle volume, in. (m ) 

3 3 

fuel-injector manifold volume, in. (m ) 

q 

single- orifice manifold or dual-orifice secondary manifold volume, in. 

(m 3 ) 

ratio of secondary (axial) to total oxidizer flow 

single- orifice pump discharge or dual- orifice pump discharge plus primary 
manifold volume, in. (m ) 

angular frequency, rad/sec 

fuel flow rate, lbm/sec (kg/sec) 

full thrust fuel flow rate, lbm/sec (kg/sec) 

oxidizer flow rate, lbm/sec (kg/sec) 
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WOP primary (tangential) flow rate, lbm/sec (kg/sec) 

WOS secondary (axial) flow rate, lbm/sec (kg/sec) 

WOSV secondary flow specified by vortex characteristic, lbm/sec (kg/sec) 

WOl full thrust oxidizer flow rate, lbm/sec (kg/sec) 

WT input value of molecular weight of combustion products, lbm/(lb)(mole) 
(g/(g)(mole)) 

WTO oxidizer molecular weight, lbm/(lb)(mole) (g/(g)(mole)) 

WTOT total propellant flow rate, lbm/sec (kg/sec) 

WTT molecular weight coefficient (ref. 8) 

X throttling index 

XG partial derivative of chamber pressure with respect to oxidizer flow, sec 
((N)(sec)/(kg)(m 2 )) 

-2 -2 

XP primary-element pressure-drop coefficient, in. (m ) 

2 2 

XS secondary-element pressure-drop coefficient, in. (m ) 

Y substituted parameter in solution of eq. (28) 

YY variable fuel-area index 


APPENDIX E 


COMPUTER LISTINGS FOR DIGITAL CHUGGING ANALYSIS 


MAIN PROGRAM 


COMMON /ABC/ W , CUM * HE SOP »COP*ELJ»RESQJ»WO»PC 
DIMENSION GAM(2), TKI2), WT(21 

DIMENSION UPLIIOO), ACRL(IOO). UPHIIOO), ACRH(IOO), FLI1D3), FH( ID 
10 ) 

DIMENSION ACCLI10). UPPL(IO), FFL(IO), ACCH(IO), UPPHI1D), FFH( 10) 
C READ IN AND WRITE VALUES OF SPECIFIC HEAT RATIO , COMBUSTION TEMP., 

C AND MOLECULAR Ml. OF PROOUCTS AT TWO CHAMBER PRESSURES FOR NOMINAL 
C MIXTURE RATIO. FOR THIS CASE N=1 IS FOR PC=8Q0 AND N=2 IS FOR PC 
C =50 0 

READ (5,11 (GAMIN! ,TKIN) ,WT(N) ,N=1,2) 

WRITE (6,1! (GAMIN). TKIN) ,WT(N) ,N=1,2) 

1 FORMAT I6X,F6.4,1X,F5.0,1X,F6.3,1X,F6.4,1X,F5.0,IX,F6.3> 
f. MM=0 OFNOTES DUAL ORIFICE 

C MM = 1 DENOTES SINGLE ORIFICE 
M M=0 
G C= 32. 2 
U F=12. 

G = 32 . 2 

C THE FOLLOWING DATA IS FOR METHANE 
R F= 12 .* 9o . 37 

C THE FOLLOWING DATA IS FUR 82.6 PERCENT FLOX 
T CR =26 1 . 

W T0 = 3f>. 8 
H V= 78. 0 
R FO = .0 5 38 
B 0 = 7 . 8 4F *4 
S T =67 . 0 
V I S = . 2 2 

R STVC=( (.041 1/RH0>*( ST/13.2)*! VIS/.19) )**.25 
W T T— I WTO/ 100.)**. 35 
H VT= I HV/1 40 . )* * . 8 

C THF FOLLOWING COEFFICIENTS ARE SUPPLIED FOR OXIDIZER PUMP. Cl IS 
C A FLOW COEFFICIENT AND C2 IS A HEAD RISE COEFFICIENT 
C 1 = 44. 253 
C 2=3.25 3F— 6 

C THF FOLLOWING ARE FULL THRUST OPERATING CONDITIONS 
IF I MM. ED. 0) GO T3 2 
PCI = 500. 

0 F = 5 • 2 5 
WO 1=10. 31 
WF1 =W.J 1 /OF 
T F 1 = 1 3 73. 

T 0= 165. 

P S = 40 . 

GO TO 3 

2 P Cl =50 0. 

0 F = 5 . 2 5 

W 0 1 = 10 . 46 
W F 1 = WO 1 /O F 
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TF1=1344. 

T 0= 156.5 
P S = 40. 

C THE FOLLOWING IS GEOMETRY DATA. A CONICAL NOZZLE IS ASSUMED AND 
C ITS LENGTH IS BASED ON SPECIFIED CHARACTERISTIC LENGTH AND 
C CONTRACTION RATIO AND THROAT AREA 

3 IF (MM.EO.OJ GO TO 4 
A T= 5 • 309 

EC-4. 

ELSTR=37. 

6LCT=1 3. 

GO TO 5 

4 AT=5.336 
£C=4. 

ELSTR=40. 

ELCT=I3. 

5 ELN=3.* <ELSTR-ELCT*EC I /(l.-2.*£C+EC**.5l 
E LC = EL G T— ELN 

V N= AT* ELN* ( l .+EC+EC**.5) /3. 

A C= AT* EC 

VC=AC*ELC 

S=VN/AC/ELN 

£LT=ELC/fcC** .44+.83*tLN/S**.33 /EC**. 2 2 
C THE FOLLOWING IS SPECIFIED INJECTOR GEOMETRY 
IF (MM.EJ.O) GO TO 6 
D 0= . 06 94 
EL0=.694 
E NF = 6 T . 

GO TO 7 

6 DP=. 01403 
DS = . 0920 
DO=.0846 
E NF-68* 

ELU= I. 

7 F C= FNE / AC 
FNO=£NF 
tNF=ENE 

C FOR DUAL DR I F I Ct INJECTOR* VOM IS MANIFOLD VOLUME FOR SECONDARY FLOW 
C AND VOS IS SUM OF MANIFOLD VOLUME FOR PRIMARY FLOW AND VOLUME BETWEEN 
C PUMP AND INJECTOR. FOR SINGLE ORIFICE INJECTOR* VOM IS INJECTOR 
C MANIFOLD VOLUME. VOS IS VOLUME BETWEEN PUMP AND INJECTOR 
C VOF IS FUEL MANIFOLD VOLUME 
IF (MM.EU.O) GU TO 8 

V 0F=4. 

FLLS=18. 

D SU = 2 . 5 
vns=5o. 

V 0M = 2 1 . 

A OV l = . 2 29 
C DO = . 6 
C CP = .6 
C CS= .6 
C CF= .6 
S WC= .448 
A PFO = l .15 8£— 2 
GO TO 9 

B AOV l =. 2 15 

0D0 = .6 
C CP — . 6 
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r. ds=.6 

C CF = .4247 

SWC=.448 

vof=4. 

ELL S = 1 8. 

0SU=2. 5 
V0S=50. 

V0M=2l. 

APEO=l .0535E-2 

C AT THIS POINT THE DECISION IS MADE WHETHER TO GENERATE BOUNDARIES 
C AT VARIOUS CHAMBER PRESSURES, USING OPERATING POINT DELAYS OR TO 
C GENERATE BOUNDARIES USING DELAY VALUES SENSITIVE TO FUEL ANNJLJS 
C AREA. JJ=0 DENOTES FORMER, JJ=1 DENOTES LATTER 

9 J = 1 

J J=1 
M N= I 
L 1 = 0 
II =0 

IF IJJ.EU.O) MN =10 
DO 89 1=1, MN 
X = I 

PC=PCI*( 1 . l-X*. II 

K =0 

A PF = AP FO 

C SPECIFY FUEL TEMPERATURE ANO CONTROL VALVE RELATIONSHIPS TO 
C CHAMBFR PRESSURE 

TF=TF]*(PCl/PC)**.il6 

A 0 V=A I V 1*1 -.009 7 7+ 1.035* (PC /PCI) -.0252*1 PC/PCI 1**2 ) 

C COMPUTE THEORETICAL CHARACTERISTIC VELUCITY AT GIVEN CHAMBER PRESSURE 
C AND MIXTURE RATIO ALSO PARTIAL W.R.T. MIXTURE RATIO 
CALL CSTRR I PC , OF , C STAR , SLOPE J 

C AT FULL THRUST, FLOWS ARE AS SPECIFIED. AT THROTTLEO CONDITIONS 
C (J NOT FOUAL TU 1), FLOWS BASED ON PREVIOUSLY COMPUTED 
C EFFICIENCY. 

IF (J.FO.ll GO TO 10 
GO TO 11 

10 W0=W01 
W.F=W£ 

GO TO 12 

11 EFF=EFF 1 
WTOT=PC*AT*GC/CSTAR/EFF 
W 0=QF/ ( OF + 1 . )* W TO T 
wF=WO/OF 

GO TO 13 

C AT FULL THRUST. COMPUTE EFFICIENCY 

12 EFF=PC*AT*GC/CSTAk/( WO+WF) 

C INTERPOLATE OR EXTRAPOLATE LINEARLY TO DETERMINE SPECIFIC HEAT 
C RATIO, MOLECULAR WEIGHT AND COMBUSTION TEMPERATURE. USE RESULTS 
C TU CALCULATE CHAMBER PRESSURE AT INJECTOR FACE. 

13 GAMMA = GAM ( 2» F< GAMI 1 ) -GAM ( 2 ) ) * I PC-500 . » /300. 

TC=1.8* I TR I 2 )«■< TK< 1I-TKI 2) )*< PC-500. ) /300. 1 
RC=15a 4./{ WT 12 IF(WT( 1 )-wT(2) ) * (PC-500. I / 300 . I 
ACOUS=SURT( GC*GAMMA*RC*TC)*UF 

E XP = < GAMM A+ 1 . ) /I GAMMA- 1.1 /2. 

FMC = ( 1 ./( 1. + (GAMMA- 1. I /2. )>**E XP/EC 
V EL C=F ML* ACO US 
F XP 1 =GAMM A / 1 GAMMA-1. ) 

PCI=PC*( 1 .*GAMMA*EMC**2) / ( 1 • F ( GAMMA- 1. ) /2. *EMC**2 ) **EXPi 
C COMPUTE VELOCITIES AND PRESSURE ORUPS 
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14 V F=WF/ENF/PC I*RF*TF/APF 

P F=PCI*SJRT< 1. -»- VF **2 /CDF **2/RF/TF /GC/UF) 

VO=WO/ENU/RHO/ 3.1416/00**2*4. 

EMR=VF/VU/0F 

DPJ = ( W0/EN0/CDQ/3. 1416/00**2*4. )**2 /RHO/2. / GC/UF 

IF IMM.FD.O) GO TO 16 

OPS=(WO/AOV)**2/RHD/2. /GC/UF 

OPP = DP J *DP S 

PC=PCI +OPP 

GO TO 17 

15 X S=4 . /ENO /COS/3. 1416/0 S* *2 

X P = 4 . / EN3 /3. /COP/3. 1416/DP**2 
W0P=W0/(1.+XP/1 l./AOV**2*XS**2 »**.51 
WOS=WO-WOP 
FOP=WOP/WO 

P SUP=PC I«-DP J+I WOS*XS 1**2 /RHO/2 . /GC/UF 
DPS = { WOP*XP I ** 2 /RHO/2. /GC/UF 
DPP = DP J +0P S 
PO=PCt «-DPP 

V OR TC=- 2. 45* (PO/PSUP ) **2*4. 52* (P0/PSUP1-1. 12 59 
WOSV=WO*VORTC 

IF I AH SI WOSV— WO SI.GT..01) GO TO 16 
GO TO 17 

16 OS=OS*( W0SV/W0S)**.5 
GO TO 15 

f. TEE FOLLOWING SECTION COMPUTES DROP SIZE. BURNING LENGTH, 

C EFFECTIVE LENGTH. DELAY TIMES, PERCENT VAPORIZED AND 
C PERFORMANCE LOSS DUE TO INCOMPLETE VAPORIZATION. AT FULL THRJST, 

C TEE PERFORMANCE LOSS DUE TO INCOMPLETE MIXING IS INFERRED. AT 
C THROTTLED CONDITION. OVERALL EFFICIENCY IS CALCULATED AND COMPARED 
C WITH THE VALUE USED. IE ERROR EXISTS, FLOWS ARE ADJUSTED 

17 RMC TN=.5* . 2 36* R ST VC* SWC*DO* SQRT(OF*VO/VF ) 

T T = ( 1.-T0/TCR»**.4 

FL50=2.75*EC** . 44* T T * I R MC TN /. 003 ) **1 . 45 * I VO/ 1 2 00 . ) ** . 7 5 *W T T *HV T/ ( P 
1C/300. »** .66 

FLFFF=FLT/EL 50*2. 7 5*EC**. 44* WT T*H VT 

T AUV =EL 50 /VO 

EIGFF=ELEFF/WTT/HVT 

PV=PPI FLGEE I 

0 EN E =C STAR*! OF + 1. I 

CALL CSTRR I PC , P V*OF , C STA A . SLO PA I 
F NUM=C S TA A* ( PV*0F+1. ) 

FTAV=ENUM/DENE 
IF (J.GT.ll GO TO 18 
FTAM=EFF/ETAV 
GO TO 20 

18 EFF=ETAV*FTAM 

IF I A8S ( EFF 1— E FF I .GT. . 001 ) GO TO 19 
GO TO 20 

19 W TOT = W0 +w F 
WTOT=WTOT*EFFi/EFF 
WO=OF/I OF + 1. )* W TO T 
WE= WO/OF 

F FF l =EFF 
GO TO 14 

C COMPUTE GAS RESIDENCE TIME, THRUST PER ELEMENT, GAS PHASE DELAY 
C COEFFICIENT and operating point pressure DROP RATIOS 

20 T HE T AG = EF F*C ST AR*EL STR /GAMMA /RC /TC/3 86 . 

TPF=10.*PC/ENE 
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OPU=DP J/PC 
OPF— ( P F~P C I ) /PC 

COEFF=l.E+3# (ELC+ELN— EL50) /VELC 
T AUM=TM ( CQEF F) 

T AU T= T A UM + T A UV 

C THE FOLLOW ING SECTION COMPUTES LINEARIZED RESISTANCES, I NDJCF ANC ES , 

C CAPACITANCES AND GAINS. 

X G=( CSf AR +( OF* !• ) * SLOPE ) /A T /GC *£FF*PCl/PC 
F G=( CSTAR—OF*< OF + 1. ) *SLOPE ) /AT /GC *EF F*PCI /PC 
R ESF=< PF* *2- PC 1**2 ) /WF/PC I 
RFSOJ = 2.*< DPP-DPS) /WO 

CALL RE SXP (Cl ,C2*RHQ. WO,PS,PO , SPEED , RE SOP) 

R ESOS=2.*DPS/WO 

ELJ=ELO*4. /ENO / 3. 1416/00**2 /UF /G 
E LSU— FLL S* 4 • /3 . 1416/USU**2/UF/G 
C OP =RHO*Vfi S/tiO * G/GC 
COM=RHG*VCM/BO*G/GC 
CUF=VOF/RF/TF 

C AT THIS POINT, PARAMETERS OF INTEREST ARE WRITTEN OUT. 

WRITE (6,21) 

2 1 FORMAT (2X,2HPC,4X,2HED,2X,2HNE,4X,2HNOf4X»3HF/E»5X,2HD0*6X,4HAF/E 
l, 6X,?HC*. 6X, 2HL* ,5X,3HEFF ,4X , 2 HWO, 6X , 2H WF , 4X ,2 HVO , 8X , 2 HPS , 5X, 2HMC, 
26 X • 3HPC I • 3X , 6H THE TAG ) 

WRITE (6,22) PC,ED,ENE,ENO,TPE ,DG , APF , C S TA R , E LST R , EF F , W3 , W F , V 0 , PS , 
IF MC , PC I ♦ T HE T AG 

22 FORMAT ( I X . F 5. 1 , 2X , F 2. 0 , 2 X ,F3 . 0 , 2 X ,F4 . 0 , 2 X , F 5 . 1 , 2X , F 5 . 4, 2 < , E9 . 4, 2< 
I, F6. I. 3X, F*. I, 3X,F 4.3, 3X,F 6.2, 2X,F 5.2 ,2X,F6.1 ,2X ,F7. I, 2X, F 4.3, 2X,F 
26.1, 2X, IPE9.3) 

WRITE (6,23) 

2 3 FORMAT < 3 X , 2 HP Q * 6X , 2HPF , 5 X , 3HD PF , 5 X , 3HDP 0, 4X , 3 HL50 , 4X , 3HL E F , 2X , 4HL 
1GFN, 3X. 2HP V, 6X , 4HT AUT, 4X, 5HF TE MP , 3 X , 3H MO R , 4X , 2 HD P, 4X,4HTAJM, 14X , 2H 
20S.6X, 4HETAM ) 

WRITE ( 6. 24) PO, PF ,OPF ,0P0 ,EL5U, EL EFF , ELGEE , PV ,T AUT ,TF , EMR, DP, T AU^ 
1, I J S , F T A M 

24 FORMAT ( IX . F 6. I ,2X ,F6. 1 , 1 X ,F6. 4,2X ,F5.4 ,3X ,F5. 3 , IX ,F5. 1 , 2X , F4. 1 , 2X 
I, F4.3, 2X, 1PE9. 3. 2X, 0PF5.0 ,1X,= 6. 3,2X,F 5.4,2 X,1 PE9. 3, 8X , DP F 5. 4, 4X,F 
24 .3 ) 

WRITE (6,25) 

25 FORMAT ( 3X. 3HW0P , 7X , 6H WOP/ W0, 4X, 3H AO V ,6 X ,6H0XG AI N, 4X , 3 HF J E L GAI N, 3< 
l, 7HSP F EOD X , 4 X, 8H0X I NJDPR , 4X , I OH UX SPL I TL)Pk , 3 X ,9 HOX PUMPRES , 4X , 8 H OX I N 
2J IN D, 4X • 6 HOX SC A P , 3 X • 6H0XMC A P ) 

WRITE (6,26) WOP, FOP, AG V, XG , F 3 , SPE ED , Q P J , OPS, RESOP, ELJ, COP, COM 

26 FURMAT ( 2X • F 6. 3 , 4X , F 6 . 4, 4X ,F 6. 4 , 4X ,F 6 . 3 , 4 X , F6. 3 , 4X , 1 PE9 . 3 , 3X, IPE9. 
13, 3X, IPE9. 3, 4X. IPE 9. 3 , 3 X , I PE9. 3 ,2 X . 1 PE9 . 3 , I X , 1 PE9 . 3/ / ) 

EFF 1=FFF 

IF ( ( I .LT.MNJ.AND. ( JJ.EO.l) ) GO TO 88 
IF ( K.EQ.O) .AND. (JJ.EO.l) ) GU TO 79 
C AT TFIS POINT THE FREQUENCY BANDS AND INCREMENT ARE SELECTED AND 
C THf SOLUTION OF THE CHARACTERI STIC EOUATIUN IS INITIATED 
F RFQ = 7 4.9 
J N=0 
N =0 
N N = 0 

27 FREQ=FREQ+.l 

IF ( (FREQ. GT. 300.) .AND. (FREQ. LT. 400. ) ) GO TO 28 
GO TO 29 

28 IF (N.EQ.O) JN = 1 
N 1 = N 

F R FQ = 400. 
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J M = 0 

N=n 

NN=1 

GO TO 30 

29 IF { FKE0.GT.fc50. ) GO TO 50 

30 W=6.?B32*FREQ 
I 1=1 

RX = 0. 

R X 1 =0. 

RX2=0. 

RRF=0. 

R R F 1 =0 . 

R RF 2=0 • 

DP0C2=0. 

OPOC 5=0. 

0P0C8=0. 

0PFf.=0. 

0 PF 1=0 • 

DPF2=0. 

C AT THIS POINT TERMS IN THE CHARACTERISTIC EQUATION ARE EVALUATED 
C AT SPFCIFIEU FREQUENCY 

EMAGF=— I. / l»/COF*PF/PCI 
T HFTAT = TAUT*W 
T HE T AM = T A UM * W 
THFTAV=TAUV*W 
T HE T Al. = FL J* W 
ALPhAG=THETAG*W 
A 1=C0S( THETAMI 
B 1= S IN I THETAM) 

A 2=C0S ( THETAV) 

B 2= S IN ( THETAV) 

A 3=COS ( THE TAT) 

B 3=S IN ( THETAT) 

Y =R 1 +AL PHAG + AI 

C AT THIS POINT THE SECOND AND FIRST ORDER COEFFICIENTS IN THE 
r. QUA OR AT IC EQUATION FOR CRITICAL OX REAL PART ARE EVALUATEO 
A G=Y— EMAGF /FG* I l.*ALPHAG**2 ) 
BL=-Xi.*B2-2.*XG*EMAGF/FG*lA3-ALPHAG*B3) 

IF (MM.EQ.O) GD TO 31 

C THE FOIL OWING SECTION COMPUTES THE REAL AND IMAGINARY PARTS OF 
C THF IMPEDANCE LOOKING UPSTREAM FOR A SINGLE ORIFICE CONFIGURATION 
RESOS=2.*DPS/WO 

DDX=(R*RE SUP*C0P)**2+l.*W**4*fc LSU**2 *CGP**2-2 . ♦W**2*ELS J*COP 
F X = { -W*CDP*RESOP**2*-W*ELSU-W**3*COP*ELSU**2)/ ( DDX) 

E X =R FSOP / 00 X 
DX = FX 

C X=R ESOS+EX 

RX=(-W*COM* (CX**2+DX**2)+DX) /I (H*COM*CX) ** 2 + (W*C0M*DX-1. ) **2) 
AX=CX/( tW*COM*CX) 0*2+1 W*C0M*DX-1.) **2) 

AXA=AX 
BXA=BX 
GO TO 32 

31 CONTINUE 

C THE FOLLOW ING SECTION COMPUTES THE REAL AND IMAGINARY PARTS OF THE 
C IMPEDANCE LOOKING UPSTREAM FOR A DUAL ORIFICE CONFIGURATION 
ft FSO S = 2 .* DP S / WOP 

RES0V=2.*< ROS/AOV)**?/772. /RHO/WOS 
R ESUB=2.*< OPS-I WOS/AOVI**2/772 ./RHO) /WOS 
UX=W*«ESU0*RESOV*COM 
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P X=RESOV+RESOB 

UDX=(W*RE SO P*C0P )**2+l .+ W**4*t LSU**2 *CQP**2— 2 ,*W**2*ELSJ*C0P 
OX=( -H*C0P*RES0P**2+W*ELSU-W**3*C0P*ELSU**2 ) / ( ODX) 

ENX=RFSOP /ODX 

E MX I =-W*R E SUB* COM/ ( RE SOV+RE SOB ) 

FLX=RFSOV 

AX=RESnB 

BXf=-R*RESOV*COM/< RE SOV+RE SOB) 

AKX=0X*RES0S**2/I (RESGS+PX) **2+0X**2) 

A HX= ( )X*I ( ELX*EMXI)**2+1. ) +E MX I* ( E NX **2 + OX**2 ) )/<EMXl**2*< ELX+ENX ) 
l**2+( l.+I OX*EMXI )**2) ) 

A YX= ( P X*R ESO S**2+R E SOS* ( P X**2+ OX**2 ) ) / ( ( RESOS+PX) **2+OX**2 ) 

GX=< E\JX*( ( ELX*EMXI )**2+l. ) + EL X*EMX I **2* ( ENX**2 +OX **2 ) ) / ( EMX I**2*( E 

11 X+ENX )** 2 + 1 l.+UX*EMXl i**Z) 

L) X=AHX +AK X 

CX=GX+AYX 

EX={ CX*< ( AX*BX1 )**2+l. ) +A X*B XI **2* I CX**2 +DX**2 ) ) / < BXI**2*< AX+CX ) ** 

12 *( l.+DX*BXl )**2) 

FX=(DX*((AX*BXI ) ** 2+ 1 • ) +BXI*(CX**2+DX**2 > ) / ( BX I**2*( AX+CX ) **2+( 1.+ 
1 DX*BX I )** 2 ) 

A XA = EX 
R XA = FX 

3? F MAGO = THF TAL +BXA 

C AT THIS POINT THE ZEROTH ORDER COEFFICIENT IN THE QUADRATIC 
C EQUATION FOR CRITICAL OX REAL PART IS EVALUATED. THE FOLLOWING 
C SFCTION SOLVES THE QUADRATIC AND RELATES SOLUTION TO THE 
C PRESSURE DROP RATIO 

CG=EMAGU**2*Y-EMAGU*XG*A2+2.*XG*EMAGF*ENAGO/FG*(B3+ALPHAG*A3 I -EM AG 
1F*EMAG0**2/FG*( l.+ALPHAG**2 )-XG**2*EMAGF/FG 
I F ( AG.NE .0. ) GO TO 33 
R X=— CG/BG 
RPC=RX- AXA 

IF (RPC.LE.O.) GO TO 27 
0PDC2=RPC*W0 /PC/2. 

GO TO 38 

33 ALPHA=BG**2-4.*AG*CG 
IF ( Al PHA ) 27. 34,35 

34 R X=— BG/2. /AG 
RPC=RX-AXA 

IF (RPC.LE.O.) GO TO 27 
UP0C2=RPC*W0 /PC/2. 

GO TO 38 

33 RX1=(-BG+SJRT( ALPHA) )/2. /AG 

K X2 = < -BG- SQR T( ALPHA ) ) / 2 . /AG 
R PC 1 =R X 1- AXA 
R PC2=R X2-AXA 

IF ( (RPC1 .LE.O. ) .AND. (RPC2.LE. 0. ) ) GO TO 27 
IF (RPC1.LE.C. ) GO TO 36 
If (RPC2.LE.C. ) GO TO 37 
l 1 = 2 

D POC 5=RPC l* WO/PC /2 • 

DP0C8=RPC2*W0/PC/2. 

GO TO 3*3 

36 RPC=RPC2 

D POC 2=KPC * WO /PC / 2. 

RX=RX? 

GO TO 38 

37 R PC = RP C 1 

DP0C? = RPC*W0/PC/2. 


58 



C AT 
C UX 
38 


39 


AO 
A 1 

A? 


A3 

AA 

C AT 
C IN 


45 


A6 


47 


c;o TO 38 

THIS POINT THE CRITICAL FUEL REAL PART IS DETERMINED USING 
SOLUTIONS. THE RESULT IS RELATED TO THE PRESSURE DROP RATIO 
RRF = ( f G*RX*Q l-FG*EMAGO*Al-XG*E MAGF*A3 + ALPHAG*EMAGO*EMAGF-EMAGF*RX i 
1/ I RX*A L PH AG + EM AGO- XG*B3 ) 

IF 4 RRF ) 27.27,43 

RRFl=( FG*RX 1 *6 1-FG*E MAGG* A L - XG*EMAGF * A3 + ALPHAG *E MAGO*EMAGF-EMAGF*R 

ix n/( ax i* alphag*emago-xg*b3 ) 

RRF2=4 FG*RX2*B 1-F G*E MAGQ* Al - XG *E MAGF* A3 + ALPH AG *E MAGGIE MAGE -EM AGF *R 
IX 2) /< *X2* AL P FA G + EMAGQ- XG*B3 ) 

IF (RRFl.GT.C. I GO TO 40 
R RF 1 =0 • 

0 POC 5 = 0 • 

IF (RRF2.LE.C. ) GO TO 41 
GO TO 42 

I F (RRF2.GT.C. ) GO TO 42 
R R F 2 = 0 . 

0 POC 8 = 0 . 

OPF 1 = SURT (PC I**2+WF*PCI*RRFU / PC- PC l /PC 
OP F2 = SQR T ( PC I**2 + WF*PCI*RRF2) /PC-PCI/PC 
GO TO 44 

OPFC = SURT{ PC I**2+WF*PCI*RRF)/PC-PC I/PC 

IF < 4 GPOC 2.EQ.0.1 .AND. 4DPFC.E0.0. I . AND. ( 0P0C5. EQ.O .) .AND. I DPF1 .EQ. 
10 •) .AND. 4 DP OC 8. E 0*0. I .AND. IDPF2.EO.O.) I GO TO 27 
THIS POINT THE LUW AND HIGH FREQUENCY SOLUTIONS ARE SORTED 
ORCFR OF INCREASING FUEL PRESSURE DROP RATIO 
IF (NN.FQ.U GO TO 47 
I F ( DPF l.EQ.O. ) GO TO 45 
l F I DPF 2. EQ. C. > GO TO 46 
N=N + 1 

UPL4NI-0PFI 
FL(N l = FRFU 
A CRL I \l I =QPUC 5 
N =N + 1 

U PL 4 N ) = DP F 2 
FL I N )=FREQ 
ACRLCN I =D POC F 
GO TO 27 


N=N+1 

UPL I N > =DPF2 
F L ( N ) = FKE 0 
A CRL I V ) =D POC 8 
GO TO 27 
N =N + l 

UPL 4 N ) =DP F I 
F L ( N )=FREQ 
A CRL ( N ) =DPOC 5 
GG TO 27 

IF 4 DPF 1. EU .0.1 GO TO 48 
IF 4 DPF 2. EQ.C. I GO TO 49 
N =N + 1 

UPHI Nl = OP F I 
F HI N I = FRE 0 
A CR H ( N )=OPOC 5 
N =N*1 

U PH { N ) =D PF2 
F HI N I =FkFU 
ACRHIN )=DP0C8 
GO TO 27 
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48 N=N + 1 
UPH(N)=DPF2 
FWNI = FREQ 
acrhin ) =opoce 
GO TO 27 

49 N=N+l 
UPH(N)=DPF1 
FF( N )=FREQ 
ACRH(N > =OPOC 5 
GO TO 27 

50 IF (IFREQ.GT.6 50.) .A NO. (N.EU.O)) JM=1 
N 2 = N 

IF ( ( JN.EQ.l I. AND. ( JM.fcO. II I GO TO 78 
IF (JN.EQ.l) GO TO 51 
IF (JN.EQ.l) GO TO 52 
I M= 3 

GO TO 53 

51 I M= 1 

GO TO 53 
5? I M=2 

GO TO 5 fa 

53 N 1 1 =N l— l 

DO 55 I X= 1 * N 1 1 

I Xl=IX+l 

00 5 5 JX= IX 1 *N1 

IF ( UPL I IXI— UPL(JX) ) 55*55*54 

54 T EMP=iJPL ( I X ) 

0 PL < IX) =UPL ( JX ) 

UPl ( JX)=TFMP 
T FMP =ACRL 1 IX ) 

A CRL ( IX ) = A C R I ( JX) 

ACRL ( JX)=TEMP 
T FMP =FL (IX) 

FL< IX)=FL<JX) 

FL( JX ) =TFMP 

55 CONTINUE 

IF ( H.EQ. 1 ) GU TO 59 

56 N ?2 = N 2— 1 

DO 58 RX= 1 * N 22 

K X 1 =K X + l 

DO 58 LX=KX1,N2 

IF ( OPH(KX)-OPHILX) ) 58.58,57 

57 T FMP -UP HI K X ) 

U PH ( K X ) =UP H ( LX) 

0 PHI I X ) =TEMP 
T FMP =ACRH I KX ) 

ACRHUX )=ACRH( LX) 

ACRHd.X >=TFMP 
T FMP=FH( K X ) 

FH( KX ) =FH( LX ) 

FH(LX)=TFMP 

58 CONTINUE 

C AT THIS POINT, THE VARIABLE ARtA CASE(JJ=1) CALLS FOR INTERPOLATION 
C TO FIND CRITICAL OX RF AL PART AND FREOUENCY AT THE OPERATING F J EL 
C PRFSSORE DROP RATIO. THE RESULT IS STURED. THE THROTTLING CASE 
C <JJ=n) CALLS FOR WRITING UUT SOLUTIONS. 

59 IF ( (JJ.EU.l ) . A NO . ( IM.FU.1) ) GO TO 60 
IF ( (JJ.EU. 1 ) . A NO . < 1M.E0.2) ) GO TO 62 
IF I (JJ.EQ. 1 ) . A NO . I IM.E0.3) ) GO TO 64 
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GO TO 68 

60 IF < (OFF. LT. GPL ( I) ) .OR. (OPF.GT .UPL (N1 ) J » GO TO 79 
L I=L !♦ 1 

N S=N L- 1 
00 61 L =1 * N S 

IF ( OPF.GT .UPL { L + l ) I GO TO 61 

ACCL ( L I ) = ACRL( L )«-(DPF-UPL (L > ) * (ACKLIL+l > — ACk L ( L ) )/(UPL(L«-l)-JPL<L) 
1 > 

F FI (LI )=FL (L 1+ ( OPF— UPL ( L) ) * (FL ( L+ 1 l-FL ( L ) ) / ( UPL ( L«- 1I-UPL ( L » ) 

UPPL ( L I l = OPF 
GO TO 79 

61 CONTINUE 
L 1=1. I- 1 
GO TO 79 

b? IF 1 ( OPF.LT.UPHdll.OR. ( DPF . G T . UPrll N2 ) ) > GO TO 79 
I L= IL + 1 
N S=N 2— 1 

00 63 L=l.NS 

IF ( OPF.GT. UPH< L+l) ) GO TO 63 

ACCHl IL » = ACKH(L ) ♦ ( OPF-UPH ( L I » * ( ACRH ( L* 1 1 - ACRH ( L ) ) / ( J PH ( L + 1 ) - J PH( L ) 

1) 

FFH( IL ) =FH(L I *• I OPF-UPH ( L I » * (F hi ( L+ 1 1 -F H ( L I) / ( UPH ( L+ 1 I -UPH ( L ) I 
0 PP ri( U ) =0P F 
GO TO 79 

63 CONTINUE 

1 L= II — 1 
GU TO 79 

64 IF ( (OPF. LT.UPLI 1)1. OR. (OPF.GT. UPL(Nl) ) ) GO TO 66 
L I=L I *■ 1 

MS=N1-1 

00 65 L =1 ♦ N S 

IF (0PF.GT.UPL(L*11) GO TO 65 

ACCL (I I I =ACR L( L I+IDPF-UPL (L) ) * (ACRL( L«-l I -ACRL( LI ) / (U PL (L + l ) -UPL ( L) 
II 

F FL ( L I )=FL(L)+( OPF -UPL ( L I 1 * (FL ( L* 1 1 -FL ( L) I / ( UPL ( L+ 1 ) -J PL ( L ) ) 

UPPL (l I l = OPF 
GO TO 66 

65 CONTINUE 
L I=L I- l 

66 IF ( (OPF. LT.UPHdl J. OR. (OPF.GT. UPHCN2I II GO '0 79 

1 L = I L 1 
NS=N2- 1 

00 67 LJ=1,NS 

IF ( OPF.GT. UPH(LJ+1I I GO TO 67 

A CCH( IL I = ACK H( L J I + ( DPF— UP H ( L J I I * ( ACRH ( LU + l 1 -ACRH(L J) )/(JPH(LJ+l)-J 
IP H( L J I I 

FFH( IL l=FH(LJ)+( DPF— UPH ( L J) ) * ( FH ( L J+ 1 ) -FH ( L J > ) / ( UPH ( L J +1 ) - JPH( L J ) I 
II PP H( IL > = 0PF 
GO TO 79 

67 CONTINUE 
I L = IL - 1 
GO TO 79 

C AT THIS POINT WRITE OUT ORDERED SOLUTIONS 

68 IF (IN. EO. 3 1 GO TO 72 
wRITF (6.69) 

69 FORMAT ( 6X. 9F>FR EUUENC Y.9 X . 14H0ELT AP FUEL/ PC . 6X . 1 2HDELT AP OX/PC I 
IF ( IM.E3.2I GO TO 71 

WRITE (6,70) (FL(KL) ,UPL(KL) .ACRL(KL) ,KL=1 ,N1) 

70 FORMAT ( 8X . OPF 5. 1 , ?X . I P2 E 2 0. 3 I 
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GU TO 7b 

71 WRITE (6*701 ( FHC KL ) *UPH (KL) * ACRH( KL) *KL=1 »N2I 

GO TO 78 

72 WRITE ( 6* 73) 

73 FORMAT ( 6X . 9HFRE0UENC Y* 9 X . 14HQELTAP FUE L / PC , 6X , l 2HDE LT AP 3X/PC*7X* 
L9 HER EQUENC Y* 7X* 14HDELTAP E UEL/ PC * 7X * 1 2HDELT AP OX/PC) 

N D=N 1 

IF INI .GT.N? ) N0=N2 

WRITE (6* 74) (FUKL) tUPLlKL) .ACRLIKL) *FH ( KL) , UPH ( KL ) , ACRHI KL ) , KL = 1 
1, NO) 

74 FORMAT i 8X*0PE5.1*2X*1P2E20*3.9X.0PF5.1 * 2 X , 1 P2 E20 . 3) 

IE (NI-N2) 75.78*77 

75 N E=N 1+1 

WRITE i £ i 76 ) (FH(KL) ,UPH(KL) .ACRH(KL) *KL=NE,N2> 

76 EURMAT ( 64X * CP E 5„ 1 * 2 X* 1P2E2 0. 3 ) 

GO TO 78 

77 NE=N2+1 

WRITE (6*70) ( F L ( KL ) • UPL ( KL ) * ACRL ( KL) »KL=NE*N1 ) 

C AT THIS PUNT* VARIABLF AREA VALUES ARE SPECIFIEO 
78 IHJJ-EQ.O) GO TO 88 

79 J =2 
K=K + 1 
Y Y=K 

IF (YY.GT.8.) GO TO 80 
IE | YY.EO.l . I GO TO 87 
A PE=APE*Y Y/ ( YY- 1. ) 

E EE1=EE E 
GO TO 14 

C AT THIS POINT • THE VARIABLE AREA VALUES OF CRITICAL OX REAL PART 
C AND FREQUENCY AI OPERATING FUEL PRESSURE DROP RATIOS ARE 
C CONVERTED TO PROPER VARIABLES FOR ORDERING AND WRITING 

80 IE << JL.EO.OI.AND.(LI.EO.O) J GU TO 89 
IF ( IL .EU.OJ GU TO 83 

I F (L I .EQ.O) GU TO 85 

00 81 I 7 = 1 * L I 
UPL ( 17 )=UPPL (12 ) 

ACRL ( 17 ) = ACC L ( 17) 

ELI 17 ) =FFL ( 17) 

81 CONTINUE 
N 1 =L I 

1 M= 3 
J J — 0 

DO 82 17 = 1* IL 
U PHI 17 ) =UPP H (17) 

ACRHI I 7 ) = ACC H( IZ) 

EH< 17 )=EFH( 17) 

82 CONTINUE 
N 2= IL 

GO TO 53 

83 DO 84 17=1. LI 
UPt ( 1/ ) =UP PL (IZ > 

ACRL ( 17 ) = ACC L ( IZ) 

F L ( IZ >=FFL( 17) 

84 CONTINUE 
N 1 =L I 

I M= 1 
J J = 0 

GO TO 53 
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85 


DO 86 17*1. IL 
UPHl 17 )*UPPHIl Z I 
A CR HI 17 l*ACCH( IZ) 
FHI I Z 1 = FFHI I Z ) 

86 CONTINUE 
N 2=IL 

I M=2 
J .1=0 

Gfl TO 56 

87 APF*.25*APF 
GO TO 14 

88 J *2 

89 CONTINUF 
STOP 
END 


SUBPROGRAM PVELG 


C THE FOLLOWING FUNCTION FITS PRIEMS PtRCENT VAP VS. ELGEN CURVE 
FUNCTION PP (ELGEE) 

D IMFNSIUN S( 17) . VI 17) 

DATA SI 1)/. 37. S <2» /1.0/,SI3)/l. 73 7 ,S<4) 72. 55/.SI 5) /3.42Z, St 6) 74. 60 
l/.SI 7) 76.80/. SI 8)78. 00/. St 9) /8.80/.SI 10) /12.20/.SI 11 )/ 15.807, S( 12) 
27 la. 107, SI 13 ) 724. 007, SI 14) 727. 007 , S I 1 5 ) 7 32 . 007 , S I 1 6 » /35 . 0 0 7 , S I 17)7 
340.007 

DATA VI 1)7.17. VI 2) 7.26/, V (3) 7.38/, VI 4) 7. 487, VI 5)/. 56/. VI 6 > 7.647 , VI 
17)7. 747, VI 8)7. 78/.VI9) 7. 807, VI 10) 7.877 ,V(il)/.917,VI 12)7. 937, VI 13) 
27.96 7. VI 14 ) 7.97 7.VI 15) /.98/.VI 16) /.9867 , VI 17 )7.99/ 

00 l 1 = 1, 16 

IF 1 I SI I) .Lfc .ELGEE) .AND. 1ELGE6 .LT.SI 1*1) ) ) GO TO 2 

1 CONTINUE 

2 PP*VI I )+( ELGEE- St l))*tVII+l)-VtI))/tStI+l>-StI>) 

IF IPP.GT..55) PP*. 99 

R ETUKN 
FNO 


SUBPROGRAM TMC 


C THIS MAP FITS THE EMPIRICAL MIXINS TIME CURVE 
FUNCTION TM ICOEFF) 

DIMENSION S 1 6 ) . VI6) 

DATA SIl)7.57,S12)71.07,S<3)7l.57,SI4)/1.87,SI5)72.0/,S<6)72.27 
DATA VI 1)7.15/, VI 2) 7. 357 . VI 3 ) 7 . 7 7 , VI4 ) 7 1 .075 7 » VI 5 ) 7 l .47 , V I 6 ) 7 1 .97 
DO l 1 = 1,5 

IF I € SI I) .LE.COEFF). AND. ICOEFF. LT.Stl + 1))) GO TO 2 

1 CONTINUE 

2 TM=1.F-3*IV< I )+ ICOEFF- SI I ) ) * I V 1 1 + 1 )- V 1 1 ) ) 7 I S 1 1 +1 > -S I I ) ) ) 

RETURN 

END 
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SUBROUTINE CSTRR 


r. THF FOLLOWING SUBROUTINE COMPUTES C* FROM MCBRIDE DATA FOR SELECTED 
f. PROPELLANT COMBINATION WITH FUEL AT DESIGN POINT T EMPERATJRE 
SUBROUTINE CSIRR { PC T , OF .C STAA . SLOPE ) 

DIMENSION V l 2. 7 ) . CS(2), SSI2I 

DATA VI I, ll/6992./»Vll»2l/7072./»V(l»31/7153./»V(l»4)/7152./»V( 1,5 
11/7151./, VI 1,6 1/7042. /.VI 1,71/6934./ 

DATA V< 2. II /702 7. /,V(2. 21/71 11./, V 12.3 1/7195./ ,V(2, 41/ 7196./, VI 2» 5 
11/7198./, VI 2,6 >/ 7086. / , V I 2 , 71 /6973 . / 

I N = 0 F/ .5 
J =IN-7 
V =5* IN 
7 = . 1*Y 

on i 1 = 1,2 

C SI I » = V t I . J ) ♦ ( VII , J+ll-VI I.J1 1 4I0F-Z1 /.5 
SSI I) = (V< I , J*l I — VI I . J) 1 /. 5 
1 CONTINUE 

C ST AA = C SI 1 1MCSI 21 -C SI 11 1*1 PC T-5 00.1 /300. 

S10PF=SSI ll + ISSI 2) -SSI II )* IPCT-500. ) /300. 

RETURN 

END 


SUBROUTINE RESXP 


C THE FOl LOWING SUBROUTINE COMPUTES FLOX PUMP RESISTANCE 
SUBROUTINE RESXP I C l ,C 2 ,RHO , WO , PS , PO. S PE ED * RES OP 1 
DIMENSION VI 101 

DATA VI1I/.6/, VI21/.61/.VI3)/. 595/ , V 1 4 1 /. 57/ , V 1 5 1 / .53/ , VI 5 1 /. 47/ , V 
l(7)/.4/,Vt81/.31/,V(91/.21/»VCl0l/.095/ 

D FPP =P0— P S 

CHN2=DEPP /C 2 /R HO / 1728. *14 4. 

P HN = 44 8 . 8 6*C 1/RHO* WO/ 1728. 

C HP H2=C HN 2/1 PHN )**2 
DO 5 1=2, 10 
7 = I 

F =VI I 1-CHPH2*! I Z-l. I*. 021**2 
IF (F.GT.O.J GO TO 5 
IF 1 ABSIF 1 .LF. .001 1 GO TO 3 
P N=< 7-1.1*. 02 
PP=I Z-2.1*.02 

V N= V I I ) 

VP = V< I- 1 ) 

S =1 VN-VP ) / .0 2 

1 PG=PP+.5*(PN-PP 1 
VG = VP + I PG-PP I* S 
F =VG— f. HPH 2*P G**2 

IF (ABSIF I . L E. .0011 GO TO 4 
IF IF.GT.O.I GO TO 2 
P N = P G 

V N = V G 
GO TO 1 

2 P P=PG 
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'■ IHI ■■ ■ 


V P= VG 
GO TO 1 

3 CHI -CHP H2* ( < Z- l.)*.02>**2 
PHI = I 7- l. I*. OP 

S=I V< I 1-V< 1-1) ) /• 02 

GO TO 6 

4 CHI=CHPH2*PG**2 
P H I = P G 

GO TO 6 

5 CONT ItMUF 

6 SPFED=448.86*C1/KH0/ 1728.* HO/PHI 
K ESOP-— Cl *C2*SP EE0*3 • 1 1 7* S 
RETURN 

END 


PRINTOUT OF RESULTS FOR DUAL-ORIFICE CONFIGRUATION - EFFECT OF 
FUEL AREA ON STABILITY LIMITS AT FULL THRUST 


pc 
son.n 
pi) 

6 3 1 . 3 
KlIP 
l .0*>h 


GAM(l) TfC(l) WT ( 1 ) GAM(2) TK(2) WT(2) 

1. Ih'JI 45 3 7 . 19.015 I.l 6 *t 6 4444. 18.930 


68 . 


NG 
66 . 
PF OPF 

654.4 C. 257 1 
wlJP / will 
U. 1010 


F/fc 
73.5 
OPO 
. 1001 
AJ V 

0.2150 


pc ed nf 
600.0 3. 68 . 

PIJ PF 

631.3 60S. 6 

W I IP WUP 


00 

. 0846 . 

L 50 
0. 769 
UXGA1 N 
4L.589 


AF/c 

1054F-UI 
LfcF 
26.6 


C* 

7152.5 


tFF 

.965 


LGEN 

55.6 


F UfcLGAI N 
61.033 


SPEEOOX 
2.396E 06 


L* 

60.3 

r Ajr 

2.132fc-33 

JX INJDPR 
5 . DO 3 E 01 


W3 

ID .66 
FT 6 MP M3R 
1364. 3.197 


WF 

1.99 


3XSP. 1 TDPR 
7 . 492 1 01 


V 3 

508.6 
D* TAUM 
0140 6.197E-04 

J X 3 JMPRE S 
3.2256 01 


ps mc pci r -u r as 

60. o .149 5 "'5. 3 9 . 3 7 4 - - 3 4 

3 S -TAM 

. 3920 .97* 

3 X I N J I N 3 3 X S C A P 3 X M C A 3 

6.771E-T 1.44T-T 


PC tU Wfc NO F/E 

500.0 3. 68 . 68 . 73.5 

PU PF OPF OPO 

631.3 1736.7 2.4609 .1001 

W f)P wf)P / rfu AJ V 

1.056 0. 10 10 0.2 150 

PC FO N F NO F /fc 

500.0 3. AH. 68 . 73.5 

Pi) PF OPF OPO 

631.3 972.8 0.9330 .1001 

Mi)P wUP/WU AO V 

1.056 0.1010 0.2150 

PC tu nf NO F/E 

500.0 3. 68 . 08 . 73.5 

Pi) PF OPF OPO 

631.3 750.3 0.48H1 .1001 

WOP wl)P/wU AUV 

l.05o 0.1010 0.2150 

PC FU 'lb NU F/fc 

500.0 3. 68 . 68 . 73.5 

PU PF OPF OPO 

631.3 656.9 0.2571 .1001 

WOP w(]P/wn AO V 

1.056 0.1010 0.2150 


DO AF/E C* 

0846 . 2 6 3 4t -0 2 7152.5 

L50 LEF LGEN PV 
0.282 66.8 51.3 .990 

□ XL.AIN FUELGAIN 

4 C. 589 41.033 


L* fcFF rfD *F VO PS 

40. D .965 ID. 46 1.99 508.6 40. 

TAJT FTEMP M JR OP TAUM 

1.212E-D3 1344. 12.789 .0143 6.584E-04 

SPEEOOX JX I N J OPR DXSP-UDT DXPJMPRES 

2 * 396 E 04 5.0036 01 7.492E 01 3.225E 01 


. 1 


DU AF/E C* l* EFF *0 rfF VD PS 

,0846 .5268E-C2 7152.6 40. D .965 13.46 1.99 508.6 43- 

L 50 LEF LGEN PV TAJT FTEMP M JR DP T AUM 

0.465 40.4 91.5 .990 1.559E-33 1344. 6.394 .0140 6.438E-04 

UXGAIN FUELGAIN SPEEOOX JXINJOPR JXSP.ITDPR JXPUMPR2S 
4C.589 41.033 2.3 96 E 04 5.003E 01 7.492E 01 3.225E 01 


PCI 
5 n 6 . 3 
3 S 

.3323 
3X1 NJI SID 
•771E-C3 3. 


m: pci 

.149 506.1 

3 S 

. C 0? 0 
JXI NJI NO 
6 . 7 7 1 E -n 3 3. 


T 9 [ T 4 r, 

9 . * 0 4 1 - 3 V 
c r AM 
.9 7* 

JXS CAP 
4 31 - - 5 l . ■ 


3 x m ; 1 1 

* 1 c - p ' 


- 3 4 


UU AF/E C* 

0846 . 79 0 It —02 7152.5 

1 50 LEF LGEN PV 
0.624 30.1 68.2 .990 

OXGAIN FUELGAIN SPEEOOX 
4 C. 5 8 9 61.033 2.396E 06 


L* EFF A 0 A F VO 

43.3 .965 13.46 1.99 508.6 

rAJT FTEMP M3R DP TAUM 
1.859E-03 1344. 4.263 .0140 6.312E 

OX I N J DPR JXSPTTDPR 
5.003E 01 7.492E 01 


PS 

40.3 


MC 

. 149 


00 AF/E C* 

0846 .1053E-01 7152.5 

L50 LEF LGEN PV 
0.769 24.4 55.4 .990 

UXGAIN FUELGAIN SPEEOOX 
40.589 41.033 2.396E 04 


34 

JXP JMPRE S 
3.225E 01 

PS 


L* EFF -10 A F V 3 

40.3 .965 13.46 1.99 508.6 40. 

TAJT FTEMP MOR DP TAJM 
2.132E-33 1344. 3.197 .0143 6.1S7E-04 

□XINJDPR JXSP.ITDPR 3XPJMPRES 

5.303E 01 7.492E 01 3.225E 01 


PCI 
5 35 .1 
3 S 

. 39? 0 
□ XI NJ I NO 
6.771E-03 1. 

MC PCI 

.149 506.1 

3S 

.0923 
OXI NJI N3 
6.771E-03 3. 


T -t E I AG 
9 . 10 4- 
T A M 
.9 71 

3XSCAP 3X*CAP 

4 1 l r - 3 5 l . 4 4 l - - T 


T H E T A 3. 

9 . 3 9 4 1-34 
E T A M 
.971 
JXSCAP 
4 31 E-OS 1 . 


NO F/E 
66 . 73.5 

OPF OPO 

C . 1686 .1001 

/wn ADV 


DO AF/E C* 

,0866 .1317E-G1 7152.5 

L 50 LEF LuEN PV 
0.904 20.6 47.1 .990 

UXGAIN FUELGAIN SPE 


L* tFF *0 -IF 

40.3 .9a5 13.46 1.99 

TAUT FTEMP MJR DP 
2.387E-33 1344. 2.558 .0140 


□XINJDPR 


3XSP. I TDPR 


VJ PS 

508.6 40.3 

TAlM 

6.090E-04 
0XPUMPR2 S 


MC PCI 

.149 536.1 

3 S 

. 392 3 
3X1 NJI N3 


T MtTAG 
9. *942-34 
E T AM 
.9 71 
3XSCAP 
4312-05 1. 

T HET AG 
9. 3942 -34 
TAM 
.9 73 
3X5CAP 


3 X MC A P 
-4 1 : -T 


3<m; \ n 
4 4 1 f - * 5 


l .056 


0 . 

10 10 

0.2150 

4C.589 

41.033 2 .396 E 04 

5.303E 01 

7.492E 01 

3. 225E 

01 

6.7 71E-33 3.431 2 -P > 1.441l-'5 

PC 

to 

NE 

NO 

F/E 

oo 

AF/E C* 

L* 

fcFF A 0 

-IF V3 


PS 

MC PCI T HfcTAG 

500 . 0 

3. 

68 

. 6 F. 

73.5 

.0846 

1 580F-0 l 7152.5 

40.3 

.965 13. 

46 1.99 508.6 

40.3 

.149 506.3 9.3942-34 

PU 


PF 

UPF 

OPO 

L50 

LEF LGEN PV 

TAJT 

FTEMP 

M DR DP 

TAUM 


JS c TAM 

631.3 


577 . 

1 0.1415 

.1001 

1.032 

18.2 41.3 .990 

2.628E- 

33 1344. 

2.131 .0140 

5.989E- 

04 

.0920 .973 

MlJP 


M)P /mii 

A J V 

0XGA1 N 

FUELGAIN SPEEOOX 

JX INJDPR 

DXSP. I TOP* 

DXP JMPR 

2S 

3XINJIN3 DX S CAP DXMCAP 

1 .056 


0 . 

1010 

U.2150 

40.589 

41.033 2.396E 04 

5.303E 01 

7.492E 01 

3.225E 

31 

6.771E-03 3,4*12-35 L. 4417-35 

PC 

EO 

NF 

NL 

F/E 

00 

AF/fc C* 

L* 

EFF A 3 

A F V3 


PS 

MC PCI T HFTAG 

500.0 

3. 

68 

. 6 8 . 

73. 5 

.0846 . 

1844E-01 7152.5 

40.3 

.963 13. 

48 2.30 509 

.7 

40. 3 

.149 506.3 9. 3742-34 

PU 


PF 

I3PF 

OPU 

L 50 

LtF LGEN PV 

TAJT 

FTEMP 

MOR DP 

TAJM 


3 S t T A M 

631.8 


559. 

4 0. 106 1 

. 1005 

1.156 

16.3 36.8 .987 

2.857E- 

33 1344. 

1.827 .0143 

5.891E- 

04 

.3920 .973 

»np 


WUP /*□ 

A3 V 

OXGAI N 

FUELGAIN SPEEUOX 

JX INJDPR 

JXSP. I TDPR 

OXPUMPR 

2 S 

D XI NJI N3 OXSCAP 3X M C A 3 

1 .0 58 


0 . 

1010 

0.2150 

4C.502 

40.946 2.397E 04 

5.024E 01 

7.524E 01 

3.227E 

01 

6.771E-03 3.4312-35 1.441E-35 
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PC 

EU NF NO 

F/t 00 

AF/E C* L* 

6FF 4U 

rfF VD 

PS 

HZ PCI THFTA 

son .0 

3. 6 H • 6b. 

73.6 . 0846 .21 0 76-0 1 7152.6 40.3 

.959 13. 

52 2.00 51 1 

1.7 40. 

,0 .149 505.1 9.1355 

pn 

Pr iJPF 

OPQ L 50 

LCF LGEN PV F Aj f FTEMP 

MOR OP 

r a j*i 

3S 5 r A M 

t>3? -H 

547.7 0.C828 

.1013 1.277 

14.7 13.3 .983 3.075E 

-33 1344. 

1.599 .0140 

5.7956-34 

.3920 .573 

MtlP 

HUP / W< J 

A3V OXOAiN 

FUbLbAI N SPEtDOX 

OX INJOPR 

DXSP. I TQp R 

OX* JMPR5 S 

0XINJIN3 JXSCAP 

1 -06? 

0. 10 10 

0.2160 40.33b 

40.780 2.4026 04 

5.0636 01 

7.583E 01 

3.2346 01 

6. 771E-03 3.4315-35 
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APPENDIX F 


PROGRAM USERS GUIDE 
Existing Program 

A general description of the digital stability program was given in the text. The 
purpose of this appendix is to describe, in more detail, the structure and operation of 
the program. This appendix, together with the symbol list in appendix D and the com- 
puter listings in appendix E, should serve as a guide in the use of the program. It 
should be noted that the existing program uses only U. S. customary units. The use of 
the gravitational conversion factor g is intended to simplify the conversion of the pro- 
gram to the International System (SI) of units. 

The program consists of a MAIN program, a curve-fit subprogram PVELG for com- 
puting the fraction of oxidizer vaporized within a generalized length (ref. 8), a curve-fit 
subprogram TMC for computing the gas-phase mixing delay (see fig. 3), a subroutine 
CSTRR for computing the characteristic velocity and its partial derivative with respect 
to mixture ratio, and a subroutine RESXP for computing the required pump speed and 
pump resistance to satisfy the operating values of pump flow rate and pressure rise. 

The MAIN program, in conjunction with the other subprograms and subroutines, 
carries out the required steady-state and stability limit calculations. Figure 14 shows 
a flow chart of the MAIN program. Statement numbers on the flow chart refer to corre- 
sponding statements found in the listing (appendix E). 

For the propellants and operating conditions considered in the example, the sensi- 
tivity of the specific heat ratio, molecular weight, and combustion temperature in the 
chamber to the mixture ratio is small. For this reason, values at the nominal mixture 
ratio are read in as data for two values of chamber pressure (500 and 800 psia) using a 
format described in statement 1. Linear interpolation is used in the MAIN program to 
compute these variables at throttled conditions. 

For the assumption of a gaseous fuel, supplied from a choked turbine and bypass 
valve, only the fuel injection temperature and gas constant are required as input fuel 
properties. For the oxidizer, the density, bulk, modulus, viscosity, critical tempera- 
ture, etc. , are required. 

Depending on the configuration to be studied, a set of operating conditions (full 
thrust chamber pressure, mixture ratio, flows) are read into the program. For the 
example, an index MM was used to identify the type of injector and associated conditions. 

The next step in the program is the selection of the type of study to be performed, 
that is, either throttling with a particular injector (JJ = 0) or varying the injector geom- 
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Figure 14. - Flow chart for main program. 









































etry (fuel area in the example) at fixed thrust conditions (JJ = 1). 

For varying injector geometry, the existing program must throttle down from full 
thrust to the desired thrust level before calculating the stability limits. 

For each specified chamber pressure at the desired mixture ratio, the subroutine 
CSTRR is called upon to compute the theoretical c*. For the full-thrust nominal- 
geometry case, the combustion efficiency is computed from the specified chamber pres- 
sure and flows, together with the calculated c* . For off-design conditions, a previously 
calculated efficiency is used to compute the flow rates. 

Using the calculated values of specific heat ratio, temperature, and molecular 
weight, the chamber pressure at the injector face is computed from the specified value 
at the throat. The injector face pressure allows one to calculate the injection velocities, 
vaporization time delay, and performance loss due to incomplete vaporization (state- 
ments 14 and 17). For the full-thrust, nominal geometry case, g performance loss due 
to incomplete mixing is inferred from the specified overall efficiency and the calculated 
vaporization loss. For off-design conditions, the mixing loss is assumed equal to the 
full-thrust value and the overall efficiency is computed and compared with the previously 
assumed value. If an error exists, the flows are readjusted and the velocities, vapori- 
zation time delay, and performance loss due to incomplete vaporization, are recomputed. 
A suitable mixing model, if available, could be inserted into the pregram to compute the 
mixing loss and overall efficiency for all operating conditions. 

System pressure drops are then computed. For the dual-orifice case, the second- 
ary (axial) restriction size is adjusted at reduced thrust levels to give a flow split that 
satisfies a specified vortex characteristic. For the example studied, the characteristic 
used was fit by the following equation: 


w 


sec 


w. 



+ 4. 52 

P sec 


1. 1259 


(FI) 


Linearized resistances, inductances, and capacitances are then calculated at each 
operating point (statement 20). The subroutine RESXP is called to compute the re- 
quired pump speed and pump resistance. The pump characteristic (^ as functioning 

r 

<pp) must be supplied by the user. For the example, values of (CHI) were read in 
for a range of <p (PHI), from 0 to 0. 18 in steps of 0.02. The subroutine RESXP uses 
linear interpolation to determine the speed required to satisfy the operating values of 
flow, discharge pressure, and inlet pressure. At this point, the operating point values 
are written out (see appendix E). 

For the case of throttling at constant fuel injector area (JJ = 0), the next step is 
the assumption of a chugging frequency or range of frequencies to be investigated (state- 
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ment 27). Experience has shown that a lower frequency band extending from about 
3. 5/r Hz (where r is the sum of the oxidizer vaporization time, gas-phase mixing de- 
lay and the gas-residence time in ms. ) to about 4. 0/ T hertz and a higher frequency band 
extending from about 8. 7/r to 12. 0/ T hertz will generally enclose the solutions of inter- 
est. Depending on the configuration and operating conditions, solutions in either or both 
bands may be found. For the examples studied, the solutions occurred at frequencies 
quite close together. For this reason, a frequency increment of 0. 1 hertz was selected 
in both frequency bands. 

At the selected frequency, the coefficients of the second-order equation (28) are cal- 
culated. First, the first- and second-order coefficients are determined. Then, the 
computer evaluates the real and imaginary parts of the oxidizer feed system impedance 
(appendix C). Using the imaginary part, the zero-order coefficient is computed (state- 
ment 32). Equation (28) is then solved for the critical values of the real part of the oxi- 
dizer feed-system impedance. If the solution(s) result in a negative value for the injec- 
tion orifice resistance, frequency is incremented and the calculations are repeated. For 
positive, real values of the oxidizer injection- orifice resistance a’, the corresponding 
values of the fuel-injector resistance are calculating using equation (29). Again, 
only positive, real values are of interest. 

At each frequency where acceptable values of a’ and are found, the resistances 
are converted to critical pressure drop to chamber pressure ratios (AP^ q /P c )’ and 
(AP-j/P c )' using equations (57) and (58). These values, together with frequency, are 
stored (statement 44). After the desired range(s) of frequency are covered, the stored 
sets of values are sorted in order of increasing values of (AP^/P c )' (statement 147). 

This choice is based on experience with these configurations where boundaries may have 
double values of (APif/P c )' for a selected value of (AP io /P c )\ The ordering of points 
makes it possible to use existing computer plot routines and also makes it easier to read 
the computer output (see appendix E). At this point, the sets of points are written out 
(statement 68). The procedure is then repeated at the lower thrust levels. 

A similar procedure is followed for the fixed-thrust, variable fuel area case. How- 
ever, the dynamic portion of the program is not used until the system is throttled down 
from full-thrust conditions to the desired P c level. The dynamic portion of the pro- 
gram is then run at the desired P c and the specified mixture ratio. For the nominal 
fuel area, a set of boundary points is obtained in either the low-frequency range (state- 
ment 60), the high-frequency range (statement 62), or both (statement 64). Since the so- 
lutions are obtained using the delay values corresponding to the operating value of 
AP^/P c , interpolation is used to determine both the frequency and (AP io /P c )' corre- 
sponding to ( AP if/ p c )' = AP if/ p c - This set of points is then stored. The fuel area is 
then changed and the procedure is repeated. After the desired range of areas has been 
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covered, the resulting sets of points are sorted in order of increasing (AP^/P c )' (see 
statement 80). The results are then written out (statement 68). 


Modification of Existing Program 

The digital program, as described previously, is set up for a gaseous methane, 
liquid-flox propellant combination operating in an expander type cycle. The program 
is, however, capable of generating stability limits for any propellant combination or 
engine configuration desired. The following section is intended as a guide in changing 
the existing program to serve one’s individual needs. 

Changes in the propellant combination can be conveniently handled if the required 
input data are available. For other gaseous-fuel liquid-oxidizer combinations, changes 
in the fuel gas constant, fuel-injection temperature schedule, oxidizer properties, and 
combustion properties are required. The case of a liquid- liquid propellant combination 
can be handled if the following points are considered: 

(1) All input data now supplied for the oxidizer are supplied for both propellants. 

(2) Two pump characteristics are supplied with their corresponding subroutines for 
speed and resistance calculations. 

(3) Both feed systems are represented by impedance which are made up of resis- 
tances, inductances, and capacitances. This representation will probably result in a 
more complicated form for the real and imaginary parts of the fuel feed- system impe- 
dance. Transformation of the critical values of the real part to pressure drop ratios 
will involve subtraction of frequency- sensitive terms as is now done with the oxidizer. 

(4) Drop sizes must be computed for both propellants. Constant drop sizes, based 
on injector geometry (ref. 8), or other suitable correlations must be used. 

(5) Computation of effective lengths, vaporization time delays, percent vaporized, 
etc. , must be made for both propellants; this will require double useage of the gener- 
alized length-percent vaporized curve. 

(6) Computation of the vaporization efficiency (ETAV) must be made using equa- 
tion (52) in reference 8. In this report, equation (53) in reference 8 is used for the case 
of a completely vaporized fuel. 

(7) The solution of the characteristic equation and the generation of stability limits 
can be accomplished by redefining the following terms: THETAM = W * (TAUM 

+ TAUVO), THETAV = W * (TAUVO - TAUVF), and THETAT = W* (TAUM + TAUVO) 
where TAUVO and TAUVF are the oxidizer and fuel vaporization times, respectively. 
TAUM is the gas-phase delay as previously defined. The vaporization delays must be 
computed separately. 
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Different feed- system configurations (valves, orifices, lines, etc.) require changes 
in the section used to compute the real and imaginary parts of the feed-system impe- 
dances. In addition, the steady-state pressure-drop calculations must be adjusted. For 
the impedance calculations, the techniques described in appendix C should be used to set 
up the reduction equations in the program. 
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